%fig7EpsBsweep.m %Code used to generate the maximum height data in figure 7 of Shemilt et %al. (2022). Calls the function longwave_evolution.m at each point on the %grid of \epsilon/B space, which solves the long-wave evolution equation. %Written by James D Shemilt, 02/2022 A = 0.25; %Initial amplitude L = sqrt(2)*pi; %domain length Tend = 1000; %Integrate to this fixed end time N = 200; %Spatial Gridpoints for simulations dEps = 0.0025; %Grid sizing in the A direction dB = 0.01; %Grid sizing in the B direction %Initial and final points for A and B Eps0 = 0.1; B0 = 0.0; Epsf = 0.2; Bf = 0.2; %Form vectors of initial A-B grid Ev = Eps0:dEps:Epsf; Bv = B0:dB:Bf; %Vectors lengths lBv = length(Bv); lEv = length(Ev); %Initialise output vectors Epsout = zeros(lBv*lEv,1); Bout = Epsout; Hmaxout = zeros(lBv*lEv,1); tpout = zeros(lBv*lEv,1); i = 1; %counter %Fill in vectors of A and B for Eps = Ev Epsout(1+(i-1)*lBv:i*lBv) = Eps+zeros(lBv,1); Bout(1+(i-1)*lBv:i*lBv) = Bv; i = i+1; end %Solve evolution equation at each point on the grid parfor i = 1:lEv*lBv Eps = Epsout(i); %Set Eps B = Bout(i); %Set B tic [H,Psi,pz,tauw,tp] = longwave_evolution(A,B,Eps,1e-6,Tend,N,L,1e-13,0.7);%solve toc Hmaxout(i) = max(H(:,end)); %Store max value of solution tpout(i) = tp; %store end-time of simulation end %Save Aout, Bout, Hmaxout, tpout solution vectors. save('HmaxoutLWEpsB.mat','Hmaxout') save('BoutLWEpsB.mat','Bout') save('EpsoutLWEpsB.mat','Epsout') save('tpoutLWEpsB.mat','tpout')