To solve coupled odes and algebraic equations

My problem is as follows:

I have initial conditions for a set of odes and algebraic equations.I want matlab to evaluate odes for intervals dt and update the algebraic equations after that.This is because all these are coupled.I am attaching a rough code of it here(syntaxes may be wrong,I am complete beginner).Some equations are yet to be defined.

syms V CA t Tu mu P Tb Qu Qb me mbw tw mbl Ab Au Rf Af Vf Vfr zw za zp xb mb r V0 f Tob r_cyl Ru Abh Auh

V(t)= (V0/2) * (r-1) * (sin(CA) + (sin(CA)*cos(CA))/sqrt((1/(f*f))-(sin(CA)*sin(CA))))

zp=V/(pi*(r_cyl^2))

P=du*Ru*Tu

Vf=V-((m-me)/du)

za = min ( zp , Rf )

if Rf>r_cyl

zw=sqrt(Rf^2-r_cyl^2);

else

zw=0;

end

Vfr = pi *( r_cyl^2 * zw + ( za -zw )* Rf^2 - ( za^3 - zw^3 )/3 ) - Vb

Af = 2*pi*Rf *( za - zw )

Vb=Vb+Vfr

Ab=Af

Rb=Rf

diff(me,t) = Af * du * ( ut * (1-exp(-(CA*60/(2*pi*rpm))/Tob)) + Sl )

me(0)=0

diff(mb,t) = Af * du * Sl + (me - mb) / Tob

diff(mu,t) = -1*diff(mb,t)

mb(0)=0

Abh = ( (pi*B^2)/2 + 4*V/B ) * (mb/( mb + mu )^0.5

Auh = ( (pi*B^2)/2 + 4*V/B ) * ( 1 - (mb/( mb + mu )^0.5 )

diff(Qu,t) = Abh * 2.43* (vp^(⅓)) *( P * Tu)^0.5 *(Tu - Twg )

diff(Qb,t) = Auh * 2.43* (vp^(⅓)) *( P * Tb)^0.5 *(Tb - Twg )

diff(Tu,t) = (diff(Qu,t)+(((mu*Ru*Tu)/P)*diff(P,t))) / (mu*cpu)

diff(Tb,t) = (diff(mb,t)*(hu-hb) + diff(Qb,t) + (((mb*Rb*Tb)/P)*diff(P,t))) / (mb*cpb)

diff(Vu,t) = diff(V,t) - diff(Vb,t)

diff(P,t) = (P((diff(mb,t)/db)+ (diff(mu,t)/du)-diff(V,t)) + ((diff(Qu,t)*Ru)/cpu) + ((diff(Qb,t)+(diff(mb,t)*(hu-hb)))*(Rb/cpb))) / ( V - ((Vu*Ru)/cpu) - ((Vb*Rb)/cpb))

Sl = Sl0 * ((Tu/Tref)^alpha) * ((P/Pref)^beta) * (1-2.06*yr^0.77)

alpha = 2.4 - 0.271*phi^3.51

beta = -0.357+0.14*phi^2.77

Tob = lt / Sl

ut = 0.08 * Ui * (du/di)^0.5

Where Ui = integral from eqn 2.43

lt = 0.8*Lv* (di/du)^0.75

diff(mbl,t) = mbw * exp(-((CA*60/(2*pi*rpm))-tw) / Tob)

Differentiating by Temp

cv = dU/dT

cp = dH/dT

Use integral from 2.70 to find h and U

gas

cp_bar/R = ai1 + ai2*T + ai3*T^2 + ai4*T^3 + ai5*T^4

h_bar/(R_bar * T) = ai1 + (ai2 / 2)*T + (ai3 / 3)*T^2 + (ai4 / 4)*T^3 + (ai5 / 5)*T^4 + ai6 / T

Fuel

Cp_bar = af1 + af2 * t + af3*t^2 + af4*t^3 + af5/(t^2)

h_bar = af1 * t + (af2 /2) * t^2 + af3 * t^3/3 + af4 * t^4 /4 - af5/t + af6 + af8

Unburned

O2

cp_bar/R = .36256*10 + -.18782/100 *Tu + .70555/(10^5) *Tu^2 + -.67635/(10^8) *Tu^3 + .21556/(10^11) *Tu^4

h_bar/(R_bar * Tu) = .36256*10 + ((-.18782/100) / 2)*Tu + ((.70555/(10^5)) / 3)*Tu^2 + ((-.67635/(10^8)) / 4)*Tu^3 + ((.21556/(10^11)) / 5)*Tu^4 + (-0.10475*10000) / Tu

N2

cp_bar/R = 0.36748*10 + -0.12082/100 *Tu +0.23240/(10^5)*Tu^2 + -0.63218/(10^9)*Tu^3 + -0.22577/(10^12)*Tu^4

h_bar/(R_bar * Tu) = 0.36748*10 + ((-0.12082/100) / 2)*Tu + (0.23240/(10^5)) / 3)*Tu^2 + ((-0.63218/(10^9)) / 4)*Tu^3 + ((-0.22577/(10^12)) / 5)*Tu^4 + (-0.10612*10000) / Tu

Fuel Kcal/g mol

Cp_bar = -0.55313+ 181.62 * ( Tu /1000 ) + -97.787*( Tu /1000 )^2 + 20.402*( Tu /1000 )^3 + -0.03095/(( Tu /1000 )^2)

h_bar = -0.55313 * ( Tu /1000 ) + (181.62 /2) * ( Tu /1000 )^2 + -97.787* ( Tu /1000 )^3/3 + 20.402* ( Tu /1000 )^4 /4 - -0.03095 / ( Tu /1000 ) + -60.751 + 20.232

Burned

O2

cp_bar/R = .36220*10 + -.73618/1000 *T - .19652/(10^6) * Tb ^2 + 0.26202/(10^10) * Tb ^3 - .28946/(10^14) * Tb ^4

h_bar/(R_bar * Tb ) = .36220*10 + ((-.73618/1000) / 2)* Tb + (- .19652/(10^6)) / 3)* Tb ^2 + ((0.26202/(10^10)) / 4)* Tb ^3 + ((- .28946/(10^14)) / 5)* Tb ^4 + (-0.12020*10000) / Tb

N2

cp_bar/R = .28963*10 + .15155/100 * Tb - .57235/(10^6) * Tb ^2 + .99807/(10^10) * Tb ^3 - 0.65224/(10^14) * Tb ^4

h_bar/(R_bar * Tb ) = .28963*10 + ((.15155/100) / 2)* Tb + (- .57235/(10^6)) / 3)* Tb ^2 + (( .99807/(10^10)) / 4)* Tb ^3 + (( - 0.65224/(10^14)) / 5)* Tb ^4 + (-0.90586*1000) / Tb

Fuel

Cp_bar = af1 + af2 * ( Tb /1000 ) + af3*( Tb /1000 )^2 + af4*( Tb /1000 )^3 + af5/(( Tb /1000 )^2)

h_bar = af1 * ( Tb /1000 ) + (af2 /2) * ( Tb /1000 )^2 + af3 * ( Tb /1000 )^3/3 + af4 * ( Tb /1000 )^4 /4 - af5/( Tb /1000 ) + af6 + af8