Assignment3 key to Introduction to MATLAB(MIT OPEN COURSEWARE)

hh.m

function r=hh()
[T,Y]=ode45(@rigid,[0 0.2],[0.5 0.5 0.5 -60]);
figure;
plot(T,Y(:,4));
ySS=Y(end,:);
for v=1:10
hold on;
[T,Y]=ode45(@rigid,[0 0.2],[ySS(:,1:3),ySS(:,4)+i]);
plot(T,Y(:,4),'r');
end

end
function dy=rigid(t,y)
C=1;
Gk=36;
Gna=120;
Gl=0.3;
Ek=-72;
Ena=55;
El=-49.4;
dy=zeros(4,1);
dy(1)=(1-y(1))*alphan(y(4))-y(1)*betan(y(4));
dy(2)=(1-y(2))*alpham(y(4))-y(2)*betam(y(4));
dy(3)=(1-y(3))*alphah(y(4))-y(3)*betah(y(4));
dy(4)=(-1/C)*(Gk*(y(1)^4)*(y(4)-Ek)+Gna*(y(2)^3)*y(3)*(y(4)-Ena)+Gl*(y(4)-El));
end

ir.m

function r=hh()
[T,Y]=ode45(@rigid,[0 0.2],[0.5 0.5 0.5 -60]);
figure;
plot(T,Y(:,4));
ySS=Y(end,:);
for v=1:10
hold on;
[T,Y]=ode45(@rigid,[0 0.2],[ySS(:,1:3),ySS(:,4)+i]);
plot(T,Y(:,4),'r');
end

end
function dy=rigid(t,y)
C=1;
Gk=36;
Gna=120;
Gl=0.3;
Ek=-72;
Ena=55;
El=-49.4;
dy=zeros(4,1);
dy(1)=(1-y(1))*alphan(y(4))-y(1)*betan(y(4));
dy(2)=(1-y(2))*alpham(y(4))-y(2)*betam(y(4));
dy(3)=(1-y(3))*alphah(y(4))-y(3)*betah(y(4));
dy(4)=(-1/C)*(Gk*(y(1)^4)*(y(4)-Ek)+Gna*(y(2)^3)*y(3)*(y(4)-Ena)+Gl*(y(4)-El));
end

escapeVelocity.m

function n=escapeVelocity(z0,c,N)
z=z0;
for n=1:N
    z=z.^2+c;
    if(abs(z)>2)
        break;
end
end

julia.m

function M=julia(zMax,c,N)
[X,Y]=meshgrid(linspace(-zMax,zMax,500),linspace(-zMax,zMax,500));
Z=X.+(Y.*i);
M=zeros(size(Z));
for i=1:size(Z,1)
    for j=1:size(Z,2)
        M(i,j)=escapeVelocity(Z(i,j),c,N);
    end
end
imagesc(atan(0.1*M));
end