How can I input two known vectors into the solve function to solve for two unknown variables for each pair of the known vectors.

I have two known vectors (Vgs and Vds). I would like to input these vectors into solve or vpasolve. I'm hoping to solve for two unknown variables (Vf and Vmos) for each individual pair of Vgs and Vds.

VgsVector = [0.6 0.625 0.65];
VdsVector = 0 : 0.01 : 1;
Vmosout = ones(numel(VgsVector),numel(VdsVector)) .* VgsVector';
Vfout = zeros(numel(VgsVector),numel(VdsVector));

syms Vf Vmos Vgs Vds;
eq1 = Vf == Vds + Vmos^2;
eq2 = Vgs = Vf + Vmos;

for i = 1 : numel(VgsVector)
    eq4 = subs(eq2, Vgs, VgsVector(i));
    for j = 1 : numel(VdsVector)
        eq3 = subs(eq1, Vds, VdsVector(j));
        S = vpasolve( eq3, eq4, [-1]);
        Vfout(i,j) = S.Vf(1);
        Vmosout(i,j) = S.Vmos(1);
    end
end

Here is how I'm currently solving this problem, with two for loops going through each value of Vgs and Vds. I'm looking for Vmos to be an array that is Vgs rows long and Vds columns wide, and for each index it is solved for each pair from Vds and Vgs.

I've converted Equation 1 to something much simpler, but in reality its far more complicated with many layers of equations. The number of unknowns remains the same though.

I'd prefer not to do this way with the for loops as it is very slow. I've had much success converting for loops into matrices, but I'm not sure if that's possible here. If there are an optimizations or changes that would help, I'd appreciate it. I'm also not entirely sure what the best practice would be if it encounters two or more possible solutions.