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

(* Version: 1.1 *)

(* This package is free software: you can resdistribut it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   GPL version 2 is available at 
   http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
   the latest GPL license is available at
   http://www.fsf.org/licensing/licenses/gpl.html
 *)

BeginPackage["UMMathDiffEq`"]

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,...
and plots the solution for the x in the range
    {x, x0, x1}
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->Thick]"

DirectionFieldPlot::usage="DirectionFieldPlot[expression, {x,x0,x1}, 
    {y,y0,y1}, plotargs->vals]
plots the direction field for the differential equation
    y'[t] == equation
for the range of x and y values given.
  Examples:
	DirectionFieldPlot[10 - 0.1 y,{x,0,10},{y,0,10}]
	DirectionFieldPlot[1 - Sin[t w],{t,0,5},{w,-1,2}]
Note that we are plotting the direction field as a graph over
two variables, not the solution of the equation, so we always
write y, not y[x], in the expression provided to DirectionFieldPlot."

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`"]

<<VectorFieldPlots`

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

DirectionFieldPlot[equation_,xvals_,yvals_,args___]:= Module[
  {aratio, scale},
  scale = Sqrt[xvals[[3]]^2 + yvals[[3]]^2]*.035;
  aratio = If[yvals[[3]]-yvals[[2]] != 0, 
	      (yvals[[3]]-yvals[[2]])/(xvals[[3]]-xvals[[2]]), 1];
  If[ aratio > 5 || aratio < 0.2, aratio = 1 ];
  VectorFieldPlot[{1, equation}, xvals, yvals, Axes->True,
		  ScaleFactor->None, ScaleFunction->(scale&),
		  AspectRatio->aratio, 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[]
(*


^*)