euler method for solving system of ODE's 1st order

Hello. I want to solve a system of first-order equations using the Euler method.

I have written this code but it does not give me the desired answer.Please check it and tell me the problem. thanks

a=0; %initial volume
b=10; %final volume
h = 0.1; % step
N = (b-a)./h;

nsteps = (b-a)/h + 1;  % this is the number of elements in t(a:h:b) (It is = N+1)
f=zeros(1,nsteps);
g=zeros(1,nsteps);
v=a:h:b;

f = 10;              % initial condition
g = 0;               
p = 0; 
q = 10;

F = @(v,f,g,p,q) -0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
G = @(v,f,g,p,q) 0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
P = @(v,f,g,p,q) -0.2*200*(p/q)+0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
Q = @(v,f,g,p,q) f+g+p;

for i=1:N
    f(i+1) = f(i) + h*F(v(i), f(i), g(i), p(i), q(i));
    g(i+1) = g(i) + h*G(v(i), f(i), g(i), p(i), q(i));
    p(i+1) = p(i) + h*P(v(i), f(i), p(i), p(i), q(i));
    q(i+1) = q(i) + h*Q(v(i), f(i), q(i), p(i), q(i));
     
    v(i+1)=a+i*h;
end
plot(v,f,'-',v,g,'--',v,p,'-o')
%plot(v,f,v,g);

[Jan]

You forgot to mention, what the desired answer is. All we see is the posted code, but this does not allow the readers to guess, what you want to do instead.

[James Tursa]

What is the differential equation you are solving? Can you post an image of it?

[Behzad Rahmani]

How to display results in a matrix??

[Jan]

You still did not mention, why you assume that your code has a problem.

Answer 1

[Alan Stevens]

Simple Euler inaccurate for large step size. Reduce step size as in following and see if you get the output you expect:

a=0; %initial volume
b=10; %final volume
h = 0.005; % step  Need a small step-size for simple Euler!!
N = (b-a)/h;
f = zeros(1,N+1);   f(1) = 10;
g = zeros(1,N+1);
p = zeros(1,N+1);
q = zeros(1,N+1);   q(1) = 10;
v=a:h:b;

F = @(f,g,p,q) -0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
G = @(f,g,p,q) 0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
P = @(f,g,p,q) -0.2*200*(p/q)+0.7*((200*(f/q))-((200^2)/50)*(g/q)*(p/q));
Q = @(f,g,p) f+g+p;

for i=1:N
    f(i+1) = f(i) + h*F(f(i), g(i), p(i), q(i));
    g(i+1) = g(i) + h*G(f(i), g(i), p(i), q(i));
    p(i+1) = p(i) + h*P(f(i), p(i), p(i), q(i));
    q(i+1) = q(i) + h*Q(f(i), q(i), p(i));
end
plot(v,f,'-',v,g,'--',v,p,'-o')
xlabel('v'), ylabel('f,g,p')
legend('f','g','p')

This results in: