How to plot two functions together using odeplot?

Question

I'm trying to compare (graphically) approximations obtained using Euler's method and RK4-method.

ode := {diff(y(x), x) = 2*cos(x)*y(x), y(0) = 1};    
p := dsolve(ode, y(x), numeric, method = classical[rk4], stepsize = .25);     
f := dsolve(ode, y(x), numeric, method = classical[foreuler], stepsize = .25)

But can't plot them together.

plots:-odeplot([p, f], x = 0 .. 10);

Error, (in plots/odeplot) input is not a valid dsolve/numeric solution

How can I solve this problem? Besides, I want to use style=point for "euler" curve, and plot exact solution on the same graph, if it's possible. What is the best method to obtain exact solution?

Any help would be appreciated.

Answer 1

The key to answering you question is that you can produce the plots separately (using odeplot or plot), and then combine them together using the plots:-display command.

For fun, let's plot the two numeric methods with style=point, at the x-values that match the fixed step-size of the forward-Euler method.

I find it more convenient to use plot instead of plots:-odeplot here, after using the output=listprocedure option of dsolve (numeric).

restart;

ode := {diff(y(x), x) = 2*cos(x)*y(x), y(0) = 1}:

stpsz := .25;

                     stpsz := 0.25

a,b := 0, 10;

                     a, b := 0, 10

numpts := floor((b-a)/stpsz + 1);

                      numpts := 41

p := dsolve(ode, y(x), numeric, method = classical[rk4],
            stepsize = stpsz, output=listprocedure):

f := dsolve(ode, y(x), numeric, method = classical[foreuler],
            stepsize = stpsz, output=listprocedure):

e := dsolve(ode, y(x)):

Pe := plot(eval(y(x),e), x = a..b, style=line, legend="Exact"):

Pp := plot(eval(y(x),p), a..b, color=blue,
           style=point, symbol=circle, symbolsize=10,
           adaptive=false, numpoints=numpts, legend="RK4"):

Pf := plot(eval(y(x),f), a..b, color=green,
           style=point, symbol=diagonalcross, symbolsize=10,
           adaptive=false, numpoints=numpts, legend="For.Euler"):

plots:-display(Pe, Pp, Pf);

You can issue the commands eval(y(x),e) and eval(y(x),f) separately, to see that they are a syntax for picking off the RHS expression or procedure from the dsolve solution.