The solution vector is assumed to be generated from computations so there is no functionality to explicitly set it during the solution process. You can prescribe an initial solution vector to the solvers with the 'init' argument.

The solution vector format is briefly explaned here "

https://featool.com/doc/feastruct":

"After solving a problem the sol field will contain the solution column vector in sol.u with rows corresponding to the degrees of freedom (ordered according to the order of the dependent variables in fea.dvar and fea.eqn.dofm). For time dependent and eigenvalue problems the columns in u correspond to solutions at different times/eigenvalues, and additionally the output times will be stored as a vector in the sol.t or sol.l field, respectively."

Basically, the rows in "sol.u" consists of the degrees of freedom for the dependent variables in blocks corresponding to the order in "fea.dvar". In the simple case of using 1st order conforming shape/basis functions ("fea.sfun = sflag1"), the location of the degrees of freedom (DOF) will correspond to the mesh vertices ("fea.grid.p' "), so that the coordinates for the degrees of freedom (corresponding to the rows in sol.u) would be equivalent to "repmat(fea.grid.p',length(fea.dvar),1)". For higher order discretizations the degrees of freedom also include values on edges, faces, and cell centers to it becomes more complicated, easiest is to use "evalinit" to get a corresponding coordinate mapping "x0" for the DOF locations (in sol.u):

n_dvar = length(fea.dvar);
x0 = [];
for i=1:length(fea.sdim)
x0 = [x0, evalinit(fea, repmat(fea.sdim(i), 1, n_dvar))];
end

Moreover, the number of degrees of freedom for each dependent variable can be found in "fea.eqn.ndof", and the actual DOF to grid cell mapping/numbering in "fea.eqn.dofmap". From the same section:

"dofm is an array specifying the local to global degree of freedom numbering for each dependent variable (size n_ldof, n_c). The rows correspond to local degrees of freedoms on each cell and the columns give the cell numbers. (For linear conforming shape functions the dof mapping will be identical to the grid.c field.) This field is created when parseprob calls the "mapdofbdr" function.

ndof is simply an help array for the numbers of degrees of freedom for each dependent variable (equals to max(dofm(:)) and also automatically generated by mapdofbdr)."