Method # 1 - Euler (a custom-defined procedure for our uses)
DO NOT CHANGE ANYTHING IN THE FOLLOWING INPUT REGION - YOU HAVE BEEN WARNED !!!
> Euler := proc( MODEL, IC, DOMAIN, N )
> local a, b, dt, dy, f, i, tt, yy, LISTpts, vars;
> vars := (op(2,lhs(MODEL) ), op([1,0],lhs(MODEL) ) );
> f := unapply( rhs(MODEL), vars );
> a := op(1,rhs(DOMAIN)); b := op(2,rhs(DOMAIN));
> dt := (b-a)/N;
> tt := a;
> yy := rhs(IC);
> LISTpts := [ tt, yy ];
> for i from 1 to N do
> dy := evalf( f(tt,yy) ) * dt;
> yy := yy + dy;
> tt := tt + dt;
> LISTpts := LISTpts, [ tt, yy ];
> od; RETURN( [ LISTpts ] );
> end:
DO NOT CHANGE ANYTHING IN THE PRECEDING INPUT REGION - YOU HAVE BEEN WARNED !!!
Let's recall what we are working with.
> MODEL;
> IC;
> VAR;
> DOMAIN;
The number of steps to be used in Euler's method is
> Nsteps := %?;
>
The points generated by Euler's method are
> ptsE := Euler( MODEL, IC, DOMAIN, Nsteps );
>
These data points can be plotted using the command
> plotE := plot( ptsE, color = GREEN, thickness = 2 ):
> plotE;
>
Another way to view the approximate solution is to overlay the computed solution on the slope field..
> display( [ plotSLOPE, plotE ] );
>