fzero() is designed for single functions of one variables that return scalar values.
fsolve() from the Optimization toolbox can handle multiple variables (and multiple functions.)
Sometimes what you can get away with is
fminunc( @(tyl) eqdiff(tyl(1), tyl(2), tyl(3)).^2, initial_values)
function dydt = eqdiff(t,y,lambda)
That last line confuses me. You have lambda on input but you assign 1 to it inside the function? What are you expecting that would do for you?
Your function has trivial solutions: just let y or lambda be 0.