Problem with Forward Euler function?

I have created a function using the forward Euler's method to approximate the solution to the ODE dy/dx(x)=f(x,y(x)). The only problem is that the function only works when y initial is a single value. I need it to work for when y initial is an array of values, but I'm not sure how I'm supposed to fix it. Here's my code:

function [x,y] = forwardEuler(f,a,b,n,y0)
h = (b-a)/n;
x = [a zeros(1,n)];
y = [y0 zeros(length(y0),n)];
for i = 1:n
    x(i+1) = x(i)+h;
    y(i+1) = y(i)+h*f(x(i),y(i));
end
end

I figured that the part that needs to be modified is where y is defined in the for loop, since y(i+1) assumes y will be a single column vector, so I tried changing it so that it counts for the actual size of y, but I'm kind of lost there.

Answer 1

[Cedric]

You should build a simple case study outside of this function. It seems that you want to build an array whose rows contain the evolution of each component of y0. You should try to see (e.g. in the command window) whether you are able to initialize such an array, and then whether you are able to fill it in correctly.

To illustrate the process:

 >> a = 2 ;   b = 12 ;  n = 10 ;  h = (b-a)/n ;
 >> y0 = [3, 4, 5] ;
 >> x = [a zeros(1,n)] ;
 >> y = [y0 zeros(length(y0),n)];
 Error using horzcat
 Dimensions of matrices being concatenated are not consistent.

well, it seems that when the user defines y0 as a row vector, it doesn't work. Therefore, we can use linear indexing to guarantee that we are dealing with a column vector:

 >> y = [y0(:), zeros(length(y0),n)];
 >> size(y)
 ans =
     3    11
 >> y
 y =
     3     0     0     0     0     0     0     0     0     0     0
     4     0     0     0     0     0     0     0     0     0     0
     5     0     0     0     0     0     0     0     0     0     0

which seems to be working. Next, we define a simple f for the test:

 >> f = @(x,y) 2*x + 0*y ;

Then we try to perform the first iteration..

 >> i = 1 ;
 >> x(i+1) = x(i)+h;
 >> x
 x =
     2     3     0     0     0     0     0     0     0     0     0

this seems to be working.. let's see how it goes for y:

 >> y(i+1) = y(i)+h*f(x(i),y(i));
 >> y
 y =
     3     0     0     0     0     0     0     0     0     0     0
     7     0     0     0     0     0     0     0     0     0     0
     5     0     0     0     0     0     0     0     0     0     0

this doesn't work; the 4 from the initial condition was replaced with 7. But if we look a little closer, it seems that you are indexing y, which is a 2D array, as if it were a 1D array. Looking at the array, it seems that step i should define y(:,i+1) based on x(i) and y(:,i) .. and here I let you bring the update. Note that the upper boundary of the range defining the loop statement is n in your code, so the last loop will define the n+1-th column of x and y, and you might want hence to review this upper boundary.

I can't write the solution for you as I suspect that it is a homework, but this shows you how to proceed: simple steps where you display/test all the steps of the computation.

[Joe]

Ok I think I got it, thanks for the help!

[Cedric]

You're welcome; don't forget to [ Accept the Answer ] if I answered your question, or ask if you need more help (keeping in mind that we/I can't give you directly the solution).

[Jan]

@Cedric: Your answer is clean, clear and descriptive, as usual. Frequently, I read even answered question to check, if I can suggest an improvement or to vote the answer and/or question. But for your answers I can omit the "suggest improvements" part and I appreciate this. (Sorry for this flattery. I do not think you need sweeties, but I do criticize often enough in this forum. I hope you join the editor team as soon as possible.)

[Cedric]

@Jan: Thank you for the positive feedback! I have to admit that I spend sometimes a little too much time on answering than what would be reasonable with my schedule, but .. I guess that most of us are in the same situation here ;-) Looking forwards to be an editor actually, if only because I am maniac about code formatting!