fsolve issue with initial condition for ODE

I am trying to solve the following ODE:

function [eta, sol] = compressible_similarity_wo

    global Gamm Ma Pr omega;

    Gamm = 1.4;
    Ma = 2;
    Pr = 0.7;
    omega=0.76;

    global eta_max_ode;

    eta_max_ode = 20;

    opt = optimset('Display','off','TolFun',1E-20);
    F = fsolve(@(F) eval_boundary(F),[0,0,0.4,1,0],opt);

    [eta_ode, fg_ode] = solve_ode(F);

    sol = [fg_ode];
    eta = eta_ode;

end

function [eta_ode, fg_ode] = solve_ode(F)
    global eta_max_ode

    options = odeset('RelTol',1e-9,'AbsTol',1e-9);
    [eta_ode, fg_ode] = ode45(@BLFunc,[0,eta_max_ode],F,options); 
    end

    function [g] = eval_boundary(F)
    % Get the solution to the ODE with inital condition F
    [eta_ode, fg_ode] = solve_ode(F);
    % Get the function values (for BCs) at the starting/end points
    f_start = fg_ode(1,1); %f(0) = 0
    df_start = fg_ode(1,2); %f'(0) = 0
    df_end = fg_ode(end,2); %f'(inf) - 1 = 0
    t_end = fg_ode(end,4); %T(inf) - 1 = 0
    dt_start = fg_ode(1,5); %T'(0) = 0
    % Evaluate the boundary function
    g = [f_start
         df_start
         df_end - 1
         t_end - 1
         dt_start];
end

function [df] = BLFunc(f)
    global Gamm Ma Pr omega;

    df = zeros(5,1);
    df(1) = f(2);
    df(2) = f(3);
    df(3) = -f(1)*f(3)/(f(4)^(omega-1))-(omega-1)*f(3)/f(4);
    df(4) = f(5);
    df(5) = -Pr*f(1)*f(5)/(f(4)^(omega-1)) - Pr*(Gamm - 1.0)*Ma*Ma*f(3)*f(3) - (omega-1)*f(5)/f(4);

end

but fsolve returns the following problem

Error using BLFunc
Too many input arguments.

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 solve_ode (line 5)
[eta_ode, fg_ode] = ode45(@BLFunc,[0,eta_max_ode],F,options);

Error in eval_boundary (line 3)
[eta_ode, fg_ode] = solve_ode(F);

Error in compressible_similarity_wo>@(F)eval_boundary(F) (line 15)
F = fsolve(@(F) eval_boundary(F),[0,0,0.4,1,0],opt);

Error in fsolve (line 230)
            fuser = feval(funfcn{3},x,varargin{:});

Error in compressible_similarity_wo (line 15)
F = fsolve(@(F) eval_boundary(F),[0,0,0.4,1,0],opt);

Error in launch (line 3)
[eta, sol] = compressible_similarity_wo;

Caused by:
    Failure in initial objective function evaluation. FSOLVE cannot continue.

Do you have an idea of what's going on?

I'll cite you the friendly manual page

The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector dydt of data type single or double that corresponds to f(t,y). odefun must accept both input arguments, t and y, even if one of the arguments is not used in the function.

That is, you simply need to change to

function [df] = BLFunc(t,f)

to get a result (no guarantee that it is THE result).