MATLAB how do I plot my Runge-Kutta4 ODE solution on whole interval

Question

I made this code for ODE solutions using Runge-Kutta4:

function y=myODE(f,y0,x0,h,x_final)
x=x0;
y=y0;
while x<x_final
  h=min(h,x_final-x);
  k1=f(x,y);
  k2=f(x+h/2,y+h*k1/2);
  k3=f(x+h/2,y+h*k2/2);
  k4=f(x+h,y+h*k3);
  y=y+(h/6)*(k1+2*k2+2*k3+k4)
  x=x+h;
end

f is the function y' = f(x,y), y0 is initial value, x0 is where the function starts, h subinterval and x_final is where the function stops.

I tried my code and it solves ODEs for me correctly, but I also want to plot my function over a xy-axis on the interval x0 to x_final with h subintervals. When I try to plot it using plot(x0:h:x_final,y) I only get an empty graph. I understand (guessing) that I have to bind my y to several x in order to plot, but how can I do that without changing my code too much?

How can I plot the graph for y given y0, interval x0 to x_final and given h?

New to MATLAB, appreciate all help I can get!

Edit: To make clear what my code is for;

I need this ODE solver both for solution and graphing. I'm supposed to study the truncation error by looking at values of y on h compared to 2*h and the stability of Runge-Kutta4 by looking at graphs of y with different h.

Answer 1

This is not a very smart refactoring of your code (because it will slow down the solving, also will kill you graphics depending on how many steps you have in your ODE) but I'm sleepy so I go for the hack:

    function y=myODE(f,y0,x0,h,x_final)
            hold(axes('Parent',figure),'on');
            x=x0;
            y=y0;
            plot(x,y, 'o');
            while x<x_final
                    h=min(h,x_final-x);
                    k1=f(x,y);
                    k2=f(x+h/2,y+h*k1/2);
                    k3=f(x+h/2,y+h*k2/2);
                    k4=f(x+h,y+h*k3);
                    y=y+(h/6)*(k1+2*k2+2*k3+k4);
                    x=x+h;
                    plot(x,y,'o');
            end;
    end

Maybe tomorrow I'll re-write this answer to something proper, but—for now—good-night! :-)

Answer 2

function y=myode(f,y0,x0,h,x_final)
x=x0;
y=y0;
plot(x0,y0,'.')
hold on
while x<x_final
  h=min(h,x_final-x);
  k1=f(x,y);
  k2=f(x+h/2,y+h*k1/2);
  k3=f(x+h/2,y+h*k2/2);
  k4=f(x+h,y+h*k3);
  y=y+(h/6)*(k1+2*k2+2*k3+k4);
  x=x+h;
  plot(x,y,'.')
  disp([x,y])
end

The comment box didn't let me to put my fixed-code in "code-format" so post it as an answer.

Answer 3

I would suggest you store the x and y values to vectors and plot them outside the loop. You may want to also output bigX and bigY in order to compare with the exact solution.

function y=myODE(f,y0,x0,h,x_final)
% Example:
%     f = @(x,y) (x.^2+cos(y));
%     y_final = myODE(f,0,0,0.1,2);
x=x0;
y=y0;

bigX = x0:h:x_final;
if bigX(end)<x_final
    % Doesn't occur when x_final = n*h for some integer n
    bigX = [bigX,x_final];
end
bigY = zeros(size(bigX));
count = 1;
bigY(1) = y0;

while x<x_final
  h=min(h,x_final-x);
  k1=f(x,y);
  k2=f(x+h/2,y+h*k1/2);
  k3=f(x+h/2,y+h*k2/2);
  k4=f(x+h,y+h*k3);
  y=y+(h/6)*(k1+2*k2+2*k3+k4);
  x=x+h;
  count = count+1;
  bigY(count) = y;
end

plot(bigX,bigY,'b-o')
xlabel('x')
ylabel('y')