I am in the process of constructing a zero dimensional adiabatic parcel model during which I wish to numerically integrate the system forward in time. The system is a closed system of first order ODEs and it is my understanding that this system is stiff which limits the numerical methods that I can utilise. I am aware that there are built in solvers for stiff problems, I do not wish to use these however I will be comparing the output of my numerical solver with them.
It is also my understanding that for a stiff system of ODEs, implicit methods are required such as the BDF and Adams-moulton solvers. Whilst I have found a large selection of literature on solving non-stiff systems I have struggled to find anything that seems applicable to my problem.
My problem is an initial value problem and the initial conditions are calculated in the main script and are then placed in a state vector of the form:
y_init = [z0, P0, T0, S0, wv0, wl0, r0];
Another vector, t, for the time range also exists which is defined in the main script:
t_end = ceil(z_top/V);
t = 0:h:t_end;
y = zeros(numel(y_init),numel(t));
y(:,1) = y_init;
At present I have assumed the step size, h, is constant, however this does not have to be the case.
I have written a function which describes the system of ODEs, below is a simplified version of this function with some calculations removed (some of these calculations depend on P, T, S ect...).
function PM_dydt = PM_diff_eq(t,y,r_dry,N,kappa,V)
z = y(1), P = y(2), T = y(3), S = y(4), wv = y(5), wl = y(6), r = y(7);
dzdt = V;
dPdt = -1.0*((g*P*V)/(Rd*Tv));
drdt = (G/r)*(S-Seq);
dwldt = (4.0*pi*(rho_w/rho_ad))*N*(r^2)*drdt ;
dwvdt = -1.0*dwldt;
dTdt = -((g*V)/Cp)-(Lv/Cp)*dwvdt;
dSdt = ((P/(eps*esw))*dwvdt-((S+1.0)*(((eps*Lv)/(Rd*T^2))*dTdt)+((g*V)/(Rd*T));
PM_dydt = [dzdt; dPdt; dTdt; dwvdt; dwldt; dSdt; drdt];
You will notice that many of the derivatives depend on knowing values of another derivative, hence the system has to be solved in the order presented.
Can anyone suggest how to numerically integrate this system forward in time without using the built in solvers? Or offer any advice on how to approach this problem?
Kind regards, Liam