clc clear all close all syms J A c0 c1 tau tau_m w0 r wmax f0 sigmaJ sigmaA sigmaE Eps Alfa W=@(tau) w0*exp(r*(tau)*(1-w0/wmax)); C=@(J,tau) sigmaJ/(1+c0*tau*J); A=@(J,tau) C(J,tau)*J; CA=@(A,J,tau) sigmaA/(1+c1*(1-tau)*A(J,tau)); sE=@(tau) sigmaE*(1-tau)^Eps; eq=simplify(solve(eval(sE(tau)*CA(A,J,tau)*f0*W(tau)*C(J,tau)-1),J)) eigen=subs(diff(eval(f0*W(tau)*C(J,tau)*J*sE(tau)*CA(A,J,tau)),J),J,eq) C_m=@(J,tau,tau_m) sigmaJ/(1+c0*J*(1/2*(tau_m+tau)-1/2*(tau_m-tau)*tanh(Alfa*(tau_m-tau)))); CA_m=@(A,J,tau,tau_m) sigmaA/(1+c1*A(J,tau)*(1/2*(2-tau_m-tau)-1/2*(tau_m-tau)*tanh(Alfa*(tau_m-tau)))); fitness=eval((C_m(eq,tau,tau_m)*CA_m(A,eq,tau,tau_m)*W(tau_m)*f0*sE(tau_m))) sel_grad=simplify(eval(subs(diff(fitness,tau_m),tau_m,tau))) br=simplify(eval(subs(diff(fitness,tau_m,2),tau_m,tau))) %%%%%%%%%%%%%%