Unable to perform assignment because the left and right sides have a different number of elements; function calls

Question

I have the following Matlab simulation loop, and I keep getting the error message:

Unable to perform assignment because the left and right sides have a different number of elements.

Error in RBC (line 44) tau(n) = torqueControl(soc(n), vehSpd(n), vehAcc(n), wh, fd, gb, eng, em, batt);

GN0 = 1;
soc(1) = 0.6; %assume batt charge at start of the wltc
for n=1:length(time_s)
GN(n) = gearControl(vehSpd(n), GN0, upSpd,downSpd);
GN0 = GN(n); %asigns the current gear as per function calculation to GN0 for next shift
tau(n) = torqueControl(soc(n), vehSpd(n), vehAcc(n), wh, fd, gb, eng, em, batt);
[soc(n+1), fuelConsumption(n), unfeas(n), engPrf(n), emPrf(n), battPrf(n), vehPrf(n)] ...
    = hev_model(soc(n), [GN(n), tau(n)], [vehSpd(n), vehAcc(n)], veh,wh,fd,gb,eng,em,batt);
end

The loop should run for 1801 steps, but when GN and soc get to the 20th step, the rest of them stay in 19. I checked similar issues, but could not find a way to fix my case. Any help would be appreciated, thanks!

Below are the functions being called in the loop.

Gear conotrol:

function GN = gearControl(vehSpd, GN0, upSpd, downSpd)
% Tested schedule
if vehSpd > upSpd(GN0)
GN = GN0+1;
elseif vehSpd < downSpd(GN0)
GN = GN0-1;
else
GN = GN0;
end
end

HEV Model:

function [x_new, stageCost, unfeas, engPrf, emPrf, battPrf, vehPrf] = 
hev_model(x, u, w, veh, wh, fd, gb, eng, em, batt)
%hev_model HEV powertrain model
%
% Input arguments
% ---------------
% x : double
%   current value of the state variable(s)
% u : double
%   current value of the control variable(s)
% w : double
%   current value of the exogenous input(s)
% veh, wh, fd, gb, eng, em, batt : struct
%   data structures for the powertain components
%
% Outputs
% ---------------
% x_new : cell
%   updated value of the state variable(s)
% stageCost : double
%   stage (running) cost incurred
% unfeas : logical
%   when set to true, marks unfeasible state/control combinations
% engPrf, emPrf, battPrf, vehPrf : struct
%   data structures to return additional quantities
%
% Model details
% ---------------
% State variables
%   x(1)    Battery SOC, -
% Control variables
%   u(1)    Gear number, -
%   u(2)    Normalised engine torque
% Exogenous inputs
%   w(1)    Vehicle speed, m/s
%   w(2)    Vehicle acceleration, m/s^2

dt = 1; % s

%% Driveline (Vehicle + Final Drive + Transmission)
[shaftSpd, shaftTrq, vehPrf] = hev_drivetrain(w(1), w(2), u(1), veh, wh, fd, gb);

%% Torque Split
% Required torque (Nm)
reqTrq = shaftTrq;

engSpd = shaftSpd.*ones(size(reqTrq)); %ones(size(reqTrq)) --> gives matrix of ones size of the reqTrq
% Maximum engine torque
engMaxTrq = eng.maxTrq(engSpd);

% Torque provided by engine
engTrq  = u(2).*engMaxTrq;
% Costraints
engUnfeas = (engTrq > engMaxTrq) | (engTrq > 0 & engSpd < eng.idleSpd & u(1)~=1) | (engTrq>0 & engSpd > eng.maxSpd) | (engTrq < 0);

% Torque provided by electric motor
emTrq  = reqTrq - u(2).* engMaxTrq;

pwtUnfeas = (reqTrq<0 & emTrq>0);

%% Torque-coupling device (ideal)
emSpd = shaftSpd .* em.tcSpdRatio;
emTrq = emTrq ./ em.tcSpdRatio;

%% Engine
% Fuel and pollutants mass flow rates
fuelFlwRate = eng.fuelMap(engSpd, engTrq); % fuel flow rate, g/s
fuelFlwRate(engTrq==0) = 0;

%% EM
% Electric motor efficiency
emSpd = emSpd.*ones(size(emTrq));
emEff = (shaftSpd~=0) .* em.effMap(emSpd, emTrq) + (shaftSpd==0);
% Limit Torque
emMaxTrq = em.maxTrq(emSpd);
emMinTrq = em.minTrq(emSpd);
% Saturate regen braking torque
emTrq = max(emTrq, emMinTrq);
% Calculate electric power consumption
emElPwr =  (emTrq<0) .* emSpd.*emTrq.*emEff + (emTrq>=0) .* emSpd.*emTrq./emEff;
% Constraints 
emUnfeas = (isnan(emEff)) + (emTrq<0)  .* (emTrq < emMinTrq) +...
                   (emTrq>=0) .* (emTrq > emMaxTrq) | (emSpd > em.maxSpd & emTrq ~= 0);

%% Battery
battPwr = emElPwr;

% columbic efficiency
battColumbicEff = (battPwr>0) + (battPwr<=0) .* batt.coulombic_eff;
% Battery internal resistance
battR = batt.eqRes(x(1));
% Battery voltage
battVoltage = batt.ocv(x(1));

% Battery current
battCurr = (battVoltage-sqrt(battVoltage.^2 - 4.*battR.*battPwr))./(2.*battR);
battCurr = real(battCurr);
% Current limits
maxChrgBattCurr = batt.minCurr(x(1));
maxDisBattCurr = batt.maxCurr(x(1));
% Saturate charge current in regen braking
battCurr = max(battCurr, maxChrgBattCurr);

% New battery state of charge
x_new(1)  = - battColumbicEff .* battCurr ./ (batt.nomCap * 3600) .* dt + x(1);
% Constraints
battUnfeas = (battPwr<=0).*(battCurr < maxChrgBattCurr) + (battPwr>0).*(battCurr > maxDisBattCurr);

%% Stage cost
stageCost  = fuelFlwRate;

%% Combine unfeasibilites
unfeas = (pwtUnfeas | engUnfeas | emUnfeas | battUnfeas);

%% Powerflow (operating mode)
pwrFlw = repmat("", size(engTrq));
pwrFlw(abs(engTrq) <= 10*eps) = 'pe';
pwrFlw(engTrq > 10*eps & emTrq > 10*eps) = 'ps';
pwrFlw(engTrq > 10*eps & emTrq < -10*eps) = 'bc';
pwrFlw(engTrq > 10*eps & abs(emTrq) < engTrq.*0.01) = 'pt';

%% Pack additional outputs
engPrf.fuelFlwRate = fuelFlwRate;
engPrf.engSpd = engSpd;
engPrf.engTrq = engTrq;
emPrf.emSpd = emSpd;
emPrf.emTrq = emTrq;
emPrf.emElPwr = emElPwr;
battPrf.battSOC = x(1);
battPrf.battCurr = battCurr;
battPrf.battVolt = battVoltage;
battPrf.battRes = battR;
vehPrf.GN = u(1);
vehPrf.engTau = u(2);
vehPrf.pwrFlw = pwrFlw;
vehPrf.vehSpd = w(1);
vehPrf.vehAcc = w(2);

Torque Control:

function tau = torqueControl(soc, vehSpd, vehAcc, wh, fd, gb, eng, em, batt)
load(fullfile("Data", "transmControlData.mat"))
load(fullfile("Data", "vehData.mat"))
GN0 = 1;
GN = gearControl(vehSpd, GN0, upSpd, downSpd);
[shaftSpd, shaftTrq, vehPrf] = hev_drivetrain(vehSpd, vehAcc, GN, veh, wh, fd, gb);
if (vehPrf.reqTrq<0)
    %regen??
else
    if(vehSpd<=7) %define vpe??
        engPrf.engTrq = 0;%pure electric
    else
        if(vehPrf.reqTrq>eng.oolTrq.Values)
            engPrf.engTrq=max(eng.oolTrq.Values,vehPrf.reqTrq-em.maxTrq.Values); %power split
        else
            if(soc<0.8) %0.8??
                engPrf.engTrq=min(eng.oolTrq.Values,vehPrf.reqTrq-em.minTrq.Values); %battery charging
            else
                engPrf.engTrq=vehPrf.reqTrq; %pure thermal
            end
        end
    end
end
tau = engPrf.engTrq;
end