close all; clear all; clc

% calculate p(z) using ODE45 integrator in Matlab
% units: SI, 

z0=0;
zmax=100;
PA=1e5;
rho0=1000;
alpha=0.5;
g=9.8;

zp=linspace(z0,zmax,100);
pp0=PA+rho0*g.*zp;
pp1=PA+rho0*g.*zp+2*alpha*(rho0*g.*zp).^1.5/(3*rho0);
pp2=PA+rho0*g.*zp+2*alpha*(rho0*g.*zp).^1.5/(3*rho0)+alpha^2*g*2.*zp.*zp/6;

[z,p]=ode45(@DE,[z0 zmax],PA); 
hold on;
plot(z,p(:,1),'k-');
plot(zp,pp0,'r-');
plot(zp,pp1,'b-');
plot(zp,pp2,'g-');
xlabel('z (m)');ylabel('p (Pa)');

function [dpdz]=DE(z,p)
dpdz=zeros(1,1);
rho0=1000;
alpha=0.5;
PA=1e5;
g=9.8;
dpdz(1)=g*(rho0+alpha*(p-PA)^0.5);
end

