function [Z_select, f_T1_select,f_T2_select] = RB_linear(Z, criteria)
%%% with deterministic initial data but the pde depends on z

%% set up
global dim;

%initial data for all Z
f_0 = ID_EQ(Z);
%Prepare
[Nx,Nv,Nz] = size(f_0);
Z_select = zeros(Nz,dim);
f_T1_select = zeros(Nx,Nv,Nz);
f_T2_select = zeros(Nx,Nv,Nz);

%% first z_1
select_i =1;
% f_T2_select(:,:,1) =  VPFP_linear(f_0(:,:,select_i),phi_x(:,select_i),Z(select_i));


%% select z_n by max(residual)
for i = 1:Nz
    %record the basis and related info & removed the selected one from the original set  
    [f_T1_select(:,:,i),f_T2_select(:,:,i)] =  VPFP_linear(f_0(:,:,select_i),Z(select_i));
    Z_select(i,:) = Z(select_i,:);
    Z(select_i,:) = [];
    f_0(:,:,select_i) = [];
    
    %Get the next basis and calculate either this basis is already enough
    Y = Residual_linear(f_T1_select(:,:,1:i),f_T2_select(:,:,1:i),...
        Z_select(1:i,:),Z);
    [select_i, error_max,~] = Max_residual_linear(Y);
    
    %If previous basis haven't reaches criteria, then get one more basis
    if error_max < criteria
        break
    end 
end


%% Get the solution for these selected parameter Z_select
f_T1_select = f_T1_select(:,:,1:i);
f_T2_select = f_T2_select(:,:,1:i);
Z_select = Z_select(1:i,:);





