solving a system of ODE

so i have a pdf attached which contains all the information to solve a system of ODE.

i have written the code in a m-file as instructed in pdf :

function ydot= fossil_fuels(t,y)
p = y(1); sigmas = y(2); sigmad = y(3); alphas = y(4); alphad = y(5);
r = [0.0 0.0 1.0 4.0 5.0 8.0 10.0 10.5 10.0 8.0 3.5 2.0 0.0 0.0];
yr = [1000 1850 1950 1980 2000 2050 2080 2100 2120 2150 2225 2300 2500 5000];
f = pchip(yr,r);
d = 8.65;
u1 = 4.95 * 10^2;
u2 = 4.95 * 10^-2;
vs = 0.12;
vd = 1.23;
w = 10^-3;
k1 = 2.19 * 10^-4;
k2 = 6.12 * 10^-5;
k3 = 0.997148;
k4 = 6.79 * 10^-2;
hs = (sigmas - (sigmas^2 - k3*alphas*(2*sigmas - alphas)))^0.5/k3;
cs = (alphas - hs)/2;
ps = k4*hs^2/cs;
ydot = [(ps - p)/d + f/u1; 1/vs * ((sigmad-sigmas)*w - k1 - u2*(ps-p)/d);...
    1/vd(k1-(sigmad-sigmas)*w);(1/vs) * ((alphad - alphas) * w-k2);...
    1/vd * (k2-(alphad - alphas)*w)];
end

Can someone please check my code. I'm sort of not very confident with my defined ydot in the last part of the code.

then to solve the question i simply write the code in another m-file:

u = yr(:,1);
v = r(:,2);
plot(u,v)

Undefined function

or variable 'yr'.

Error in Untitled2

(line 1)

u = yr(:,1);

Thank you!

Answer 1

[darova]

Function looks okay except f/u1. Has to be (as instruction says):

ydot = [(ps - p)/d + ppval(f,t))/u1 ...

Read about ode45. Main .m file has to be like:

tspan = [1000 5000];                            % time span (i took from yr)
init = [1 1 1 1 1];                             % initial conditions for p, sigmas,sigmad,alphas,alphad
[t,y] = ode45(@fossil_fuels,tspan,init);

p =      y(:,1); 
sigmas = y(:,2); 
sigmad = y(:,3); 
alphas = y(:,4); 
alphad = y(:,5);