function [f_0_select, Z_select, f_T_select] = RB_init(Z, criteria)

%initial data for all
f_0 = ID_EQ(Z);
[Nx, Nv, Nz] = size(f_0); N = Nx*Nv;
f_0 = reshape(f_0, N,Nz);

%first z_1
[select_i, error_max] = Max_L2(f_0);
f_0_select = zeros(N,Nz);
Z_select = zeros(Nz,1);

for i = 1:Nz
    %record the basis and related info from the previous step & removed the selected one from the original set
    f_0_select(:,i) = f_0(:,select_i);
    f_0_select(:,1:i) = Orthogonal(f_0_select(:,1:i));
    f_0(:,select_i) = [];
    Z_select(i,:) = Z(select_i,:);
    Z(select_i,:) = [];
    %If previous basis haven't reaches criteria, then get one more basis
    if error_max < criteria
        break
    end
    %Get the next basis and calculate either this basis is already enough
    [select_i,error_max] = Max_Proj(f_0,f_0_select);
end
%solution basis
%reshape everything to matrix
f_0_select = reshape(f_0_select(:,1:i),Nx,Nv,i);
Z_select = Z_select(1:i,:);

%% This is for do projection onto the initial subspace
%Get the solution for these basis
[~,f_T_select] = VPFP_linear(f_0_select, Z_select);






