You need to initialize the arrays before loop for speed and to store the values into an array.
NOTE: Check the last values of IC,IK, gCa and gK. As of now they are zeros.
clc; clear all ;
%% Forward Euler Method
close all; clear all; clc;
%Constants set for all Methods
Cm=20; % Membrane Capcitance uF/cm^2
dt=0.01; % Time Step ms
t=0:dt:200; %Time Array ms
I=90; %External Current Applied
ECa=120; % mv Na reversal potential
EK=-84; % mv K reversal potential
El=-60; % mv Leakage reversal potential
gbarCa=4.4; % mS/cm^2 Na conductance
gbarK=8; % mS/cm^2 K conductance
gbarl=2; % mS/cm^2 Leakage conductance
% Initialize
V = zeros(size(t)) ;
minf = zeros(size(t)) ;
w = zeros(size(t)) ;
gCa = zeros(size(t)) ;
gK = zeros(size(t)) ;
ICa = zeros(size(t)) ;
IK = zeros(size(t)) ;
V(1)=-65; % Initial Membrane voltage
V1 = -1.2;
V2 = 18;
w(1)=aw(V(1))/(aw(V(1))+bw(V(1))); % Initial w-value
for i=1:length(t)-1
%Euler method to find the next m/n/h value
minf(i+1) = 1/2*(1 + tanh((V(i)-V1)/V2));
w(i+1)=w(i)+dt*((aw(V(i))*(1-w(i)))-(bw(V(i))*w(i)));
gCa(i)=gbarCa*minf(i);
gK(i)=gbarK*w(i);
gl=gbarl;
ICa(i)=gCa(i)*(V(i)-ECa);
IK(i)=gK(i)*(V(i)-EK);
Il=gl*(V(i)-El);
%Euler method to find the next voltage value
V(i+1)=V(i)+(dt)*((1/Cm)*(I-(ICa(i)+IK(i)+Il)));
end
%Store variables for graphing later
FE=V;
FEw=w;
%Plot the functions
figure
plot(t,FE);
legend('Forward Euler');
xlabel('Time (ms)');
ylabel('Voltage (mV)');
title('Voltage Change for Hodgkin-Huxley Model');
figure
plot(t,FEw,'b',t,minf,'g');
ylabel('Gaiting Variables')
xlabel('Time (ms)')
legend('w Euler','minf Euler');
figure('Name', 'Conductance')
plot(t,gCa, 'b', t, gK, 'g', 'LineWidth', 2)
legend('Action Potential', 'Ca+ Conductance', 'K+ Conductance')
xlabel('Time (ms)')
ylabel('Voltage (mV)')
title('Conduction of K+ and Ca+')
figure('Name', 'Currents')
plot(t,ICa, 'r',t,IK, 'b', 'LineWidth', 2)
legend('ICa+', 'IK+')
xlabel('Time (ms)')
ylabel('Current')
title ('Currents')
%%HH Function
%%AlphaBeta Function
function a=aw(v) %Alpha for Variable w
V3 = 2; % mV
V4 = 30; % mV
phi = 0.04; %1/ms %Rate scaling parameter
a=1/2*phi* cosh((v-V3)/(2*V4))*(1 + tanh((v-V3)/V4));
end
function b=bw(v) %Beta for variable w
V3 = 2; % mV
V4 = 30; % mV
phi = 0.04; %1/ms %Rate scaling parameter
b=1/2*phi* cosh((v-V3)/(2*V4))*(1 - tanh((v-V3)/V4));
end
