Haoyue

Bugfixes solveGenericBVP

+misc/solveGenericBVP.m

@@ -72,6 +72,24 @@ else
 		LDisc = L;
 	end
 end
+if isa(A0,'function_handle')
+	A0Disc = eval_pointwise(A0,zdisc);
+else
+	if size(A0,1)~= ndisc % constant matrix
+		A0Disc(1:ndisc,:,:) = repmat(shiftdim(A0,-1),ndisc,1,1);
+	else
+		A0Disc = A0;
+	end
+end
+if isa(H,'function_handle')
+	HDisc = eval_pointwise(H,zdisc);
+else
+	if size(H,1)~= ndisc % constant matrix
+		HDisc(1:ndisc,:,:) = repmat(shiftdim(H,-1),ndisc,1,1);
+	else
+		HDisc = H;
+	end
+end
 n = size(LDisc,2);
 T1 = [eye(n); zeros(n,n)]; % selection matrix of first solution state
 T2 = [zeros(n,n); eye(n,n)]; % selection matrix of second solution state
@@ -193,7 +211,7 @@ for jBlock = 1:length(sigMod.blockLength)
 		for zeta_idx=1:z_idx
 			temp_prod(zeta_idx,:,:) = ...
 				reshape(Psi_mu(z_idx,zeta_idx,:,:),[2*n 2*n])...
-				*T2*(reshape(LDisc(zeta_idx,:,:),[n n])\reshape(A0(zeta_idx,:,:),[n n]))*T1.';
+				*T2*(reshape(LDisc(zeta_idx,:,:),[n n])\reshape(A0Disc(zeta_idx,:,:),[n n]))*T1.';
 		end
 		A0_til(z_idx,:,:) = trapz(zdisc(1:z_idx),temp_prod,1);
 		Y(z_idx,:,:) = reshape(Psi_mu(z_idx,1,:,:),[2*n,2*n]) - reshape(A0_til(z_idx,:,:),[2*n 2*n]);
@@ -226,7 +244,7 @@ for jBlock = 1:length(sigMod.blockLength)
 		g2jk = G2*varphi_j_k;
 		
 		for z_idx = 1:ndisc
-			hjk(z_idx,:) = reshape(H(z_idx,:,:),[n nv])*varphi_j_k;
+			hjk(z_idx,:) = reshape(HDisc(z_idx,:,:),[n nv])*varphi_j_k;
 			temp_prod = zeros(z_idx,2*n,1); % preallocate and reset
 			for zeta_idx=1:z_idx
 				temp_prod(zeta_idx,:,1) = ...
@@ -323,10 +341,10 @@ for i=1:n
 	end
 end
 for z_idx=1:ndisc
-	lhs(z_idx,:,:) = sd(L(z_idx,:,:))*sd(Xzz(z_idx,:,:))...
-		+ sd(M(z_idx,:,:))*sd(sol(z_idx,:,:))...
+	lhs(z_idx,:,:) = sd(LDisc(z_idx,:,:))*sd(Xzz(z_idx,:,:))...
+		+ sd(muSys(z_idx,:,:))*sd(sol(z_idx,:,:))...
 		+ sd(sol(z_idx,:,:))*S ...
-		+ sd(A0(z_idx,:,:))*sd(sol(1,:,:));
+		+ sd(A0Disc(z_idx,:,:))*sd(sol(1,:,:));
 end
 % rhs = H;
 % err = lhs-rhs;