Undefined function or variable 'capacite_APD'

Hello everybody,

I have a problem wit matlab:

%% TP simulation de photodiode à avalanche
% version 1.0 Novembre 2018
% Louis Giraudet

%% graphics settings
    global customFontSize customLineWidth customMarkerSize
    customFontSize=14; customLineWidth=1.5; customMarkerSize=4;
    set(0,'DefaultAxesFontSize', customFontSize);
    set(0,'DefaultLineLineWidth', customLineWidth);
    set(0,'DefaultLineMarkerSize', customMarkerSize);
    global figureCounter %axisHandle thatFigure
    figureCounter=0;


%% Constantes physiques
    global h_Planck k_Boltzmann q_Elementary m_Electron permeability_vacuum light_speed epsilon_vacuum q_over_epsilon_vacuum
    h_Planck                = 6.62606957e-34;           % joule.second from NIST web site
    k_Boltzmann             = 1.3806488e-23;            % joule/kelvin from NIST web site
    q_Elementary            = 1.602176565e-19;          % coulomb from NIST web site
    m_Electron              = 9.10938291e-31;           % kilogram from NIST web site
    permeability_vacuum     = 4*pi*1e-7;                % henri/meter (exact)  
    light_speed             = 2.99792458e8;             % meters/second (définition, exact)
    epsilon_vacuum          = 1/(permeability_vacuum*light_speed*light_speed); % farad/meter
    q_over_epsilon_vacuum   = (q_Elementary/epsilon_vacuum);

%% Constantes de la structure APD
% Eclairement
    lambda = 1.55e-6; % (m) longueur d'onde d'éclairement 1.55µm
    Popt = 1e-9; % (W) puissance optique 1nW
% Courant d'obscurité
    Iobs = 1e-6; % (A) Courant d'obscurité, supposé constant
% Résistance de polarisation
    Rp = 100; % résistance de polarisation, à ajuster
% Surface de la diode
    Diametre=3e-5;
    Surface  = pi*Diametre^2/4;
% Si3N4
    n_Si3N4 = 1.98; % indice de réfraction de la couche anti-reflet
% AlInAs
    vn_aval = 4.5e4;    % (m/s) vitesse de dérive en saturation
    vp_aval = 4.5e4;    % (m/s)
    k = 0.3;            % rapport des coefficients d'ionisation
    Waval = 0.4e-6;     % (m) épaisseur de la zone d'avalanche
    epsr_AlInAs = 12.4; % permittivité relative
    n_AlInAs = 3.21;    % indice de réfraction
    tau_n_aval = Waval/vn_aval; % (s) temps de transit
    tau_p_aval = Waval/vp_aval; % (s) temps de transit
% GaInAs
    vn_abs  = 6.5e4;        % (m/s) vitesse de dérive en saturation
    vp_abs  = 4.8e4;        % (m/s)
    Wabs = 1.5e-6;          % (m) épaisseur de la zone d'absorption
    epsr_GaInAs = 13.7;     % permittivité relative
    alpha_abs = 0.8e6;      % (m-1) coefficient d'absorption à 1.55µm
    tau_n_abs = Wabs/vn_abs; % (s) temps de transit
    tau_p_abs = Wabs/vp_abs; % (s) temps de transit
    

%% Frequences et pulsation
    fmax = 100e9;                        % 100 GHz
    frequence = fmax*logspace(-3,0,301); % 3 décades en dessous de fmax
    omega = 2*pi*frequence;              % pulsation

%% Calcul de la capacité de l'APD
    % calcul de la capacité pour une valeur de Wabs
    Cj = capacite_APD(Wabs, epsr_GaInAs, Waval, epsr_AlInAs, Surface);
    fprintf('Capacité de jonction de l''APD\n\t Cj = %3.3g pF \n',Cj*1e12);
    % calcul de la capacité pour un vecteur de Wabs
    Wabs_array = linspace(1.0,2.0,6)*1e-6; % 1 à 2µm par pas de 0.2µm
    Waval = linspace(0.1,0.5,5)*1e-6; % 0.1 à 0.5µm par pas de 0.1µm
    Cj = capacite_APD(Wabs_array, epsr_GaInAs, Waval(1), epsr_AlInAs, Surface);
    legendList = build_legend_list('Waval=', Waval*1e6, 'µm');
    graphique('new', 'loglin', 'Capacité de jonction', '',  Wabs_array*1e6, 'Wabs (µm)', Cj*1e12, 'Cj (pF)', '-');
    for j=2:length(Waval)
        Cj = capacite_APD(Wabs_array, epsr_GaInAs, Waval(j), epsr_AlInAs, Surface);
        graphique('hold', 'loglin', 'Capacité de jonction', legendList,  Wabs_array*1e6, 'Wabs (µm)', Cj*1e12, 'Cj (pF)', '-');
    end
    
%% Calcul du coefficient de réflexion de l'antireflet
%% Coefficient de reflexion R:
d=linspace(0,500e-9,100);
n1=1;
A=(n1*n_AlInAs-(n_Si3N4)^2)^2*sin((2*pi/lambda)*n_Si3N4*d).^2;
B=(n_Si3N4)^2*(n_AlInAs-n1)^2*cos((2*pi/lambda)*n_Si3N4*d).^2;
C=(n1*n_AlInAs+(n_Si3N4)^2)^2*sin((2*pi/lambda)*n_Si3N4*d).^2;
D=(n_Si3N4)^2*(n_AlInAs+n1)^2*cos((2*pi/lambda)*n_Si3N4*d).^2;
R=(A+B)./(C+D);
  graphique('new', 'loglin', 'Coefficient de reflexion', '',  d, 'd(m)', R, 'R', '-');
 %% L'epaissuer optimale:
epaisseur_optimale=min(R);
%% 
r_couplage=1-min(R);
%% Sensibilité primaire Rprimaire:
%freq=light_speed/lambda;
%coefficient d'absorption du GalnAs:
Wabs=linspace(0,3e-6);
alpha_abs = 0.8e6; 
r_interne=(1-exp(-alpha_abs.*Wabs));
Rprimaire=r_couplage.*r_interne*(lambda*1e6)/(1.24);
graphique('new', 'lin', 'Rprimaire', '', Wabs, 'Wabs(m)', Rprimaire, 'Rprimaire', '-');
%% Taux de génération
Surface  = pi*(Diametre)^2/4;
G0=(Popt*Rprimaire)./(Surface*Wabs*q_Elementary);
graphique('new', 'lin', 'Taux de génération', '', Wabs, 'wabs', G0, 'G0', '-');

%% Fonction de transfert électrique
% h(w)=Ue/iph
% Ue=i0*exp(jwt)
% iph=(Rp*i0*exp(jwt)/1+j*w*Cj*Rp)
h=(Rp./(1+i.*omega.*Cj(1).*Rp));
graphique('new', 'loglin', '', ' Fonction de transfert électrique',omega, 'omega', 10*log10(abs(h)), 'h', '-');
%% Coefficients d'ionisations:
M=linspace(1,10);
k=0.3;
alpha_n=(1./((1-k)*Waval)).*log(M/(1+k.*(M-1)));
alpha_p=k.*alpha_n;
%% Les différents constantes:
taux_n_abs=Wabs/vn_abs;
taux_p_abs=Wabs/vp_abs;
a=alpha_n-(i.*omega./vn_aval);
b=alpha_p-i*omega/vp_aval;
r1=(b-a)/2-sqrt(((b-a).^2/4)+(a*b-alpha_n*alpha_p));
r2=(b-a)/2+sqrt(((b-a).^2/4)+(a*b-alpha_n*alpha_p));

%% La densité de courant total:
% Densité de courant abs:
A1=(1-exp(-i*omega*taux_n_abs))/(i*omega*taux_n_abs);
A2=1-(1-exp(-i*omega*taux_p_abs))/(i*omega*taux_p_abs);
J_abs=G0*(Wabs/(Wabs+Waval)).*(1/(i*omega*taux_n_abs).*(1-A1))+(1/(i*omega*taux_p_abs).*(1-A2))+A1.*A2.*(alpha_n*(exp(r2*Waval)-exp(r1*Waval))./((r2+a)*exp(r1*Waval)-(r1+a)*exp(r2*Waval)));
% Densité de courant aval:
J_aval=G0*(Waval/(Wabs+Waval)).*(A1).*((r2(r2+a-alpha_n)*(exp(r1*Waval)-1)-r1*(r1+a-alpha_n)*(exp(r2*Waval)-1))./(Waval*r1*r2((r2+a)*exp(r1*Waval)-(r1+a)*exp(r2*Waval))));
% Densité de courant total:
J_total=J_abs+J_aval; 
%% Fonction de transfert intrinsèque

%% Rapport Signal/Bruit
    
%%%%%%%%%%% FONCTIONS UTILES   %%%%%%%%%%%%%%%%%%%%

%% Capacité de jonction
%function [Cj] = capacite_APD(Wabs, epsr_abs, Waval, epsr_aval, Surface)
%% Calcul de la capacité de la jonction de l'APD
% comme la mise en série des deux capacités absorbant et avalanche
% auxquelles on ajoute une capacité d'effets de bord de 50fF
    %global epsilon_vacuum
    %Cbord = 50e-15; % capacité parasite sur les bords de la jonction, inévitable 50fF
   % Cabs = epsilon_vacuum*epsr_abs*Surface./Wabs;    % capacité de la couche absorbante
    %Caval = epsilon_vacuum*epsr_aval*Surface./Waval; % capacité de la couche d'avalanche
    %Cj = Cbord + Cabs.*Caval./(Cabs+Caval); % capacité totale = Cabs et Caval en série
%end

It show me an error: Undefined function or variable 'capacite_APD'.

Error in TP_APD_etudiants (line 65)

Cj = capacite_APD(Wabs, epsr_GaInAs, Waval, epsr_AlInAs, Surface);