ODE solution problem

*To*: mathgroup at smc.vnet.net*Subject*: [mg93105] ODE solution problem*From*: Virgil Stokes <vs at it.uu.se>*Date*: Sun, 26 Oct 2008 01:29:39 -0500 (EST)

I do not understand why the following code (using Mathematica vers. 5.2) does not give the correct solution. That is, the variable swt in the equation module is never set to 1, although the value of y(t) exceeds the value of yhi. The correct solution should have the following appearance, * * * * * * * * where it should reach a maximum value of 0.1 (yhi) at about 15 seconds, and then start to decrease, reaching a minimum of 0.025 (ylo) at about 30 seconds, etc. That is, the solution should oscillate between ylo and yhi. The following code however, when plotted shows a linear increase in y(t). I would appreciate it greatly if someone can point out my error(s) in the code, since I have not been able to identify why this code does not work properly. ClassicalRungeKutta /: NDSolve`InitializeMethod[ClassicalRungeKutta, __] := ClassicalRungeKutta[]; ClassicalRungeKutta[___]["Step"[f_, t_, h_, y_, yp_]]:= Block[{deltay, k1, k2, k3, k4}, k1 = yp; k2 = f[t + 1/2 h, y + 1/2 h k1]; k3 = f[t + 1/2 h, y + 1/2 h k2]; k4 = f[t + h, y + h k3]; deltay = h (1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4); {h, deltay} ]; c = 0.6; g = 9.81; R = 0.05; r = 0.007; yhi = 0.1; ylo = 0.025; Qin = 0.00005; eq[y_,swtch_] := (Qin - swtch*c*Sqrt[2*g*y]*Pi*r^2)/(Pi*R^2); Remove[equation]; swt = 0; equation[y_] := Module[{}, If[y <= ylo,swt=0;ex = eq[y,swt]]; If[y >= yhi,swt=1;ex = eq[y,swt]]; ex]; sol = NDSolve[{y'[t]== equation[y[t]],y[0]==0},y[t],{t,0,35},Method->ClassicalRungeKutta,StartingStepSize->0.02]; --V. Stokes