Choosing a numerical method for solving a stiff system of ODES

Hello,

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); % choose when to stop the solver (in this case after 250m)

t = 0:h:t_end;         % set the time range, h is the step size

y = zeros(numel(y_init),numel(t)); % initialise an empty matrix

y(:,1) = y_init;       % set the first column values to the initial conditions

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);

% 1) Define constants

% CONSTANTS ARE DEFINED HERE

% 2) Calculate parameters

% CALCULATIONS GO HERE

% 3) Define ODE system
 
  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));

  % Create a vector showing the time derivatives
  
  PM_dydt = [dzdt; dPdt; dTdt; dwvdt; dwldt; dSdt; drdt];
  
end

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