(*^
*)
(* :Title:  diffeq.m *)
(* :Author: Gavin LaRose *)
(* :Copyright: 1999 *)
(* :Summary: Differential Equations Functions *)

BeginPackage["Local`diffeq`"]

ODESolve::usage="ODESolve[{derivs==stuff, y[x0]==otherstuff,...}, y[x], {x,x0,x1}, 
plotargs->vals]  calculates the solution to the differential 
  equation 'derivs = stuff' with initial conditions 'y[x0] = 
  otherstuff,...', finding y[x] for x between x0 and x1, and plots it.  
  The number of initial conditions supplied must equal the order of the 
  differential equation.
  	Example: ODESolve[{y''[x] == -y'[x] - (y[x])^2, y[0] == 1.5, y'[0] == 0},
  y[x], {x,0,5}, PlotStyle->RGBColor[0,0,1]];"

ODEMSolve::usage="ODEMSolve[y'[x]==f[x,y], {y[x0],y0,y1,ystep}, y[x],
{x,x0,x1}, plotopt->vals] finds the solution to the first-order
differential equation y'[x] = f[x,y] for each of the initial
conditions y[x0] = y0, y[x0] = y0 + ystep,..., y[x0] = y1.  All
solutions are obtained for x0 <= x <= x1, and are plotted.  If ystep
is omitted, it is assumed to be 1.
	Example: ODEMSolve[y'[x] == x y[x] - (y[x])^2, {y[0],0,5}, y[x], {x,0,4}, 
PlotRange->{0,6}]"

Begin["`Private`"]
ODESolve[equations_,y_,xvals_,args___]:= Module[{solnrule,ysoln},
  solnrule = NDSolve[equations,y,xvals,MaxSteps->100000];
  ysoln = y /. solnrule[[1]];
  Plot[ysoln,xvals,args]
]

ODEMSolve[equation_,icvals_,y_,xvals_,args___]:= 
  Module[{icvar,icstart,icstop,icinc,numsoln,iceqn,i,eqns,solnrule,ysoln,
      allsolns},
	icvar = icvals[[1]];
	icstart = icvals[[2]];
	icstop = icvals[[3]];
 If[Length[icvals] > 3,
			icinc = icvals[[4]],
			icinc = 1];
	numsoln = Round[(icstop-icstart)/icinc]+1;
		Do[iceqn = icvar == icstart + (i-1)*icinc;
			eqns = Flatten[{equation,iceqn}];
			solnrule = NDSolve[eqns,y,xvals,MaxSteps->100000];
			ysoln[i]= y/.solnrule[[1]], {i,1,numsoln}];
		allsolns = Table[ysoln[i],{i,1,numsoln}];
		Plot[Evaluate[allsolns],xvals,args]]

End[]
EndPackage[]
(*


^*)