Hi Matlab Pros,

I need help with this code. I am using RK4 method to solve 2 non-linear equations. yprime is defined quadratically so a, b and c has to be defined. solve these ODEs store the values of y(x) and Ar(x). Equation of Theta is given, plots of x0 vs theta and Ar vs theta should be presented.

function [Theta,xo]=GM_CC
r=0.2;
alpha=50;
x_f=0.209; %inital value
x_0=0.001; %targeted value
n = 209;



Theta =zeros(n,1);
x0=zeros(n,1);

for count=1:1:n
x0(count)=count*((x0_max-x0_min)/n);
Theta(count)=GM_CC(r,x0(count),alpha);

count=count+1;
end
plot(Theta,x0,'m')



function Theta = GM_CC(r,x0,alpha)

% This function simultaneously solves equations A1.10(gives y(x)) and A1.16(gives Ar(x) using 4th order Runge
% Kutta for Counter-Current Outputs of this function include theta (cut-fraction) y, and Ar where Ar is defined as
% Am/Lo
% x0 = mol% O2 in non-permeate stream,x0=(0.209 to 0.001)



h=10^-4; % discretizes space into fine subdivisions for better approximation
xf=0.209; %standard value for oxygen content in air
x=x0+0.1*h; % At time=0, both numerator and denominator are 0, so this offset is added
r=0.2; % ratio of low side to high side pressure.
Am=5.4732; %Calculated Area of Membrane from Halvorson handout (in m^2)

% Defined variables note that the equations include yprime which is solved
% quadratically so we define a,b,c and then the yprime.
a=1-alpha;
b=-1+alpha+(1/r)+(x/r)*(alpha-1);
c=-alpha*(x/r);
yp=(-b+sqrt(b^2-4*a*c))/(2*a);
Ar=0; % Initial value
y = yp; % Stores y' and Ar as one variable
nsteps=ceil((xf-x0)/h); % Defines steps
X(1,1)=x0; Y(1,1)=y'; AR(1,1)=Ar; % Storage vectors for X,Y, and AR

% For loop to solve equations
for i=1:nsteps;
% This ensures we reach xf on last step.
if i==nsteps
h=xf-x;
end

f=rhs1(x,y,x0,r,alpha);
K1=h*f;
f1=rhs1(x+h/2,y+K1/2,x0,r,alpha);
K2=h*f1;
f2=rhs1(x+h/2,y+K2/2,x0,r,alpha);
K3=h*f2;
f3=rhs1(x+h,y+K3,x0,r,alpha);
K4=h*f3;

F=rhs2(x,y,x0,Ar,r,alpha);
J1=h*F;
F1=rhs2(x+h/2,y,x0,Ar+K1/2,r,alpha);
J2=h*F1;
F2=rhs2(x+h/2,y,x0,Ar+K2/2,r,alpha);
J3=h*F2;
F3=rhs2(x+h,y,x0,Ar+K3,r,alpha);
J4=h*F3;

y=y+(K1+2*K2+2*K3+K4)/6;
Ar=Ar+(J1+2*J2+2*J3+J4)/6;

x=x+h;

X(i+1,1)=x;
Y(i+1,1)=y;
AR(i+1,1)=Ar;
end
plot(x0,Ar(i,1),'k')
Theta = (xf-x0)/(Y(i+1,:)-x0);
%Calculates the cut fraction (theta) from final value of y
%display(Theta) %displays calculated cut fraction
%display(Lo) % displays calculated Lo value
%display(Ar) % displays calculated Ar value
%display(y) %displays calculated y value



% The following functions store values that are solved for in the for loop
% shown above
function f=rhs1(x,y,x0,r,alpha) %This function is used to solve for y using equation A1.16
if x==x0
x=x0+0.1*10^-4;
end
% Variables a, b, c, and yp change with time and must be included in this
% function.
r=0.2;
a=1-alpha;
b=-1+alpha+(1/r)+((x/r)*(alpha-1));
c=(-alpha*x)/r;
yp=(-b+sqrt(b^2-4*a*c))/(2*a);

f=((y-x0)*((1-y)*alpha*(x-r*yp)-y*((1-x)-r*(1-yp))))/((x-x0)*((1-x)*alpha*(x-r*yp)-x*((1-x)-r*(1-yp)))); %This is equation A1.16


function F=rhs2(x,y,x0,~,r,alpha) %This funciton is used to solve for Ar using equation A1.10
% Variables a, b, c, and yp change with time and must be included in this
% function.
r=0.2;
a=1-alpha;
b=-1+alpha+(1/r)+((x/r)*(alpha-1));
c=(-alpha*x)/r;
yp2=(-b+sqrt(b^2-4*a*c))/(2*a);

F=(-((y-x0)/(x-y)))/((((1-x)*alpha*(x-r*yp2)-x*((1-x)-r*(1-yp2))))); %This is equation A1.10

Please help me or give suggestions for the code.