# MATLAB: Simulating impedance using ODE45

differential equationsimaginary numbersode45

I am trying to simulate some impedance spectra. As a part of this I need to solve a system of coupled differential equations given as:
u''(x) = z1/z2 u(x) : eq(1)
and
i''(x) = z1/z2 i(x) : eq(2)
The impedance i'm looking for is in the last point: x=L
and I know that the analytical solution is:
sqrt(z1*z2)*coth(L(sqrt(z1/z2)).
My boundary conditions are:
u'(0) = 0
i(0) = 0
u(L) = A*(exp(L*sqrt(z1/z2))+exp(-L*sqrt(z1/z2)))
i(L) = A*sqrt(z2*z1)*(exp(L*sqrt(z1/z2))-exp(-L*sqrt(z1/z2)))
The code I have so far is:
``syms  C1(x) C2(x) Y % frequency resolution:FrequencyPoints = 100;% Vectro of simulated frequncies of modeling:frequencies = logspace(2.5,6,FrequencyPoints);% Preallocation of results:z_num = zeros(1,FrequencyPoints);% Run over all frequencies.for f=1:FrequencyPoints        w = frequencies(f);        % Z1 = resistor, Z2 = capacitor    Z1 = 0.1;    Z2 = 1/complex(0,w*10);        %     DC1  = diff(C1,1);    D2C1 = diff(C1,2);    DC2  = diff(C2,1);    D2C2 = diff(C2,2);        Eq1 = D2C1 == Z1/Z2*C1(x);    Eq2 = D2C2 == Z1/Z2*C2(x);        [VF,Subs] = odeToVectorField(Eq1, Eq2);    ftotal = matlabFunction(VF, 'Vars',{x,Y});    ic = zeros(2,1);    xspan = [0,1];    [x,Y] = ode45(@(x,Y) ftotal(x,Y), xspan, ic);    z_num(f) = C1(2)/C1(2);end``
The error I get when this runs is:
``Index exceeds the number of array elements (2).Error in symengine>@(x,Y)[Y(2);sqrt(1.0e+1).*Y(1).*1.0e+2i;Y(4);sqrt(1.0e+1).*Y(3).*1.0e+2i]Error in ODEexperiment>@(x,Y)ftotal(x,Y) (line 24)    [x,Y] = ode45(@(x,Y) ftotal(x,Y), xspan, ic);Error in odearguments (line 90)f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.Error in ode45 (line 115)  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);Error in ODEexperiment (line 24)    [x,Y] = ode45(@(x,Y) ftotal(x,Y), xspan, ic);``

``% frequency resolution:FrequencyPoints = 100;% Vectro of simulated frequncies of modeling:frequencies = logspace(2.5,6,FrequencyPoints);% Preallocation of results:z_num = zeros(1,FrequencyPoints);L = 1;A = 0.1;% Run over all frequencies.for f=1:FrequencyPoints    w = frequencies(f);    % Z1 = resistor, Z2 = capacitor    z1 = 0.1;    z2 = 1/complex(0,w*10);        u = A*(exp(L*sqrt(z1/z2))+exp(-L*sqrt(z1/z2)));    i = A*sqrt(z2*z1)*(exp(L*sqrt(z1/z2))-exp(-L*sqrt(z1/z2)));    z_num(f) = u/i;end``