clc
%% Question 1

%------------------------Euler Method-----------------------
dx_eu = 1e-4; % Given dx for Euler's method
beta = 0; % C(10)

Y_init =[0; %y1(0) =0, from the initial value prescribed
0; %v1(0) =0 from equations 3-4
0; % y2(0) =0
1]; %v2(0) =1

Y_eu =Y_init; %Y_eu will be updated with the values i+1 in each iteration
for x = 0:dx_eu:10

f = rhs1(x,Y_eu(:,end));
Y_ip1 = Y_eu(:,end)+f*dx_eu;
Y_eu = [Y_eu Y_ip1];

end
y1_b = Y_eu(1,end-1);
y2_b = Y_eu(3,end-1);
theta = (beta-y1_b)/y2_b;
C_eu = Y_eu(1,1:end-1)+theta*Y_eu(3,1:end-1);
figure(1)
plot(0:dx_eu:10,C_eu)
xlabel('x')
ylabel('c')
title('Solution using Euler method')

%------------------------Range Kutta Method-----------------------
dx_rk = 1e-2;
Y_rk = Y_init;
for x = 0:dx_rk:10

K = dx_rk*rhs1(x,Y_rk(:,end));
Y_ip1 = Y_rk(:,end)+dx_rk*rhs1(x+0.5*dx_rk,Y_rk(:,end)+0.5*K);
Y_rk = [Y_rk Y_ip1];

end
y1_b = Y_rk(1,end-1);
y2_b = Y_rk(3,end-1);
theta = (beta-y1_b)/y2_b;
C_rk = Y_rk(1,1:end-1)+theta*Y_rk(3,1:end-1);
figure(2)
plot(0:dx_rk:10,C_rk)
xlabel('x')
ylabel('c')
title('Solution using Range-Kutta method')

%--------------------- @ode45------------------------------
dx_ode =0.1;
X = 0:dx_ode:10;
[X,Y_ode] = ode45(@rhs1,X,Y_init);
y1_b = Y_ode(end,1);

Problem 1 Create a function M-file in MATLAB to solve the system o...