Question: Getting stuck writing my own package

Dear there,

I want to create a new package. 

Goal of the package:

I have two functions psi(n,m,x) and w(n,x).

These functions are defined differently in case A, and in case B. 

In my future works,

I want to use the functions psi(n,m,x) and w(n,x) in ONE (same) worksheet by calling the package.

Case A;

 K:=2^(k-1):
with(orthopoly): 
psi:=(n,m,x)->piecewise((n-1)/K <= x and x <= n/K,
(sqrt((m+1/2)*2^k))*LegendreP(m,2^k*x-2*n+1), 0);   
w:=(n,x)->1;  
 

Case B;

K:=2^(k-1):
with(orthopoly):
unprotect(gamma):
h:=(m,epsilon,gamma)->2^(epsilon+gamma+1)*GAMMA(epsilon+m+1)*GAMMA(gamma+m+1)/((2*m+1+epsilon+gamma)*m!*GAMMA(epsilon+gamma+1));


psi:=(n,m,x)->piecewise((n-1)/K <= x and x <= n/K,  2^(k/2) *h(m,epsilon,gamma)*simplify(JacobiP(m,epsilon,gamma,2^(k)*x-2*n+1)), 0); 
  

w:=(n,x)->(1-x)^(epsilon)*(1+x)^gamma; 

How can we achieve this in the simplest way?

MY TRY: (Please click to download the code.mw)

test:=module()
export functionA,functionB;
 option package;


######################################################
functionA:=proc(k,M) local K,h,psi,w; 
K:=2^(k-1):
with(orthopoly): 
psi:=(n,m,x)->piecewise((n-1)/K <= x and x <= n/K,
(sqrt((m+1/2)*2^k))*LegendreP(m,2^k*x-2*n+1), 0);   
w:=(n,x)->1;   
return psi(n,m,x) ,w:
end proc:
######################################################


functionB:=proc(k,M,epsilon,gamma) local K,h,psi,w; 
K:=2^(k-1):
with(orthopoly):
unprotect(gamma):
h:=(m,epsilon,gamma)->2^(epsilon+gamma+1)*GAMMA(epsilon+m+1)*GAMMA(gamma+m+1)/((2*m+1+epsilon+gamma)*m!*GAMMA(epsilon+gamma+1));


psi:=(n,m,x)->piecewise((n-1)/K <= x and x <= n/K,  2^(k/2) *h(m,epsilon,gamma)*simplify(JacobiP(m,epsilon,gamma,2^(k)*x-2*n+1)), 0); 
  

w:=(n,x)->(1-x)^(epsilon)*(1+x)^gamma;   
return psi(n,m,x),w(n,x):
end proc:
######################################################

end module;


savelib(test):


# FOR EXAMPLE
restart:
with(test);

#For k=2,M=3 in FunctionA, let's calculate values of "Psi(2)" and "w(2,x)"
k:=2:
M:=3:
K:=2^(k-1):
N:=K*M:
(psi,w):=functionA(2,3):
psi:=(n,m,x)->psi; 
w:=(n,x)->w;
unprotect(Psi):
Psi:=x->Array([seq(seq(psi(i,j,x),j=0..M-1),i=1..K)] )^+:
Psi(2);
w(2,x);


#Now, let's calculate values of "Psi(2)" and "w(2,x)" For  k=3,M=4, epsilon=1, gamma=2 in FunctionA 

 (psi,w):=functionB(3,4,1,2):
Psi(2);
w(2,x);