(*:Mathematica:: V2.0 *)

(*:Context: "ODE" *)

(*:Title: Ordinary Differential Equations *)

(*:Author: Dennis M. Snow, Department of Mathematics *)

(*:Summary: Some functions for solving ode problems in Math 226 *)

(*:History: Created Feb 20, 1992 *)

BeginPackage["Calculus`ODE`"];

DirectionField::usage=
"DirectionField[f, {x, x0, x1}, {y, y0, y1}] plots the direction field for the differential equation y' = f in the region x0 <= x <= x1, y0 <= y <= y1. The number of vectors to plot is determined by the option PlotPoints->n"

Options[DirectionField]:={PlotPoints->10};

IntegralCurves::usage=
"IntegralCurves[y'[x] = f[x,y[x]], y[x], {x,xmin,xmax}] plots integral curves of the given differential equation from xmin to xmax. Initial points can be specified with the option InitialPoints->{{x0,y0}, {x1,y1}, ...}."

Options[IntegralCurves]:={InitialPoints->{{0,-1},{0,0},{0,1}}}

Begin["`Private`"]

DirectionField[f_,{x_,x0_,x1_},{y_,y0_,y1_},opts___Rule]:=
Module[{m,n,nf},
  n=PlotPoints/.{opts}/.Options[DirectionField];
  m=Min[Abs[x1-x0],Abs[y1-y0]]/(2n);
  Show[Graphics[Table[
  nf=N[f];Line[{{x,y},{x,y}+m{1,nf}/Sqrt[1+nf^2]}],
  {x,x0,x1,(x1-x0)/n},{y,y0,y1,(y1-y0)/n}] ],
  Axes->Automatic] ]

IntegralCurves[eqn_, y_, xvals_List, opts___Rule]:=
Module[{solns, i, pts},
  pts=InitialPoints/.{opts}/.Options[IntegralCurves];
  solns = Table[NDSolve[Flatten[Join[{eqn},
    {(y/.{xvals[[1]]->pts[[i,1]]//N})==pts[[i,2]]//N}] ],
    y,xvals],{i,1,Length[pts]}];
  Plot[Evaluate[y/.solns],xvals] ]

End[]

EndPackage[]