diff --git a/matlab/swz/c-code/sbvar/switching/switch.c b/matlab/swz/c-code/sbvar/switching/switch.c
index 71db57c883acbf29136d72341a754a9780413887..008b560edec01e40611def532d7fd53950e3b754 100644
--- a/matlab/swz/c-code/sbvar/switching/switch.c
+++ b/matlab/swz/c-code/sbvar/switching/switch.c
@@ -1286,7 +1286,7 @@ int* CreateStateIndex(TMarkovStateVariable* sv, TMarkovStateVariable** list, int
   int* index;
   if (!(index=dw_CreateArray_int(sv->nstates)))
     {
-      fprintf(stdout,"CreateStateIndex():  Out of memory.\n");
+      printf("CreateStateIndex():  Out of memory.\n");
       exit(0);
     }
   for (i=sv->nstates-1; i >= 0; i--)
@@ -1294,7 +1294,7 @@ int* CreateStateIndex(TMarkovStateVariable* sv, TMarkovStateVariable** list, int
       for (k=j=0; j < n; j++)
 	if ((k=IncrementIndex(k,i,sv,list[j])) == -1)
 	  {
-	    fprintf(stdout,"CreateStateIndex():  Unable to find required state variable.\n");
+	    printf("CreateStateIndex():  Unable to find required state variable.\n");
 	    exit(0);
 	  }
      index[i]=k;
diff --git a/matlab/swz/c-code/sbvar/var/PrintDraws.c b/matlab/swz/c-code/sbvar/var/PrintDraws.c
index ac98add7586868e7d62a6cf2427838cb803df71a..a980ef0348f2476e40fbf1037533da24c7d84ccf 100644
--- a/matlab/swz/c-code/sbvar/var/PrintDraws.c
+++ b/matlab/swz/c-code/sbvar/var/PrintDraws.c
@@ -49,7 +49,7 @@ int main(int nargs, char **args)
   //=== Help Screen ===
   if (dw_FindArgument_String(nargs,args,"h") != -1)
     {
-      fprintf(stdout,"print_draws <options>\n");
+      printf("print_draws <options>\n");
       PrintHelpMessages(stdout,include_help,additional_help);
       return 0;
     }
@@ -66,7 +66,7 @@ int main(int nargs, char **args)
   dw_initialize_generator(seed);
 
   //=== Setup model and initial parameters
-  fprintf(stdout,"Reading data...\n");
+  printf("Reading data...\n");
   if (!(model=CreateTStateModelFromEstimateFinal(nargs,args,&cmd)))
     {
       fprintf(stderr,"Unable to read model or parameters\n");
@@ -89,19 +89,19 @@ int main(int nargs, char **args)
   free(filename);
 
   // Burn-in period with calibration of jumping parameters
-  fprintf(stdout,"Calibrating jumping parameters - %d draws\n",tuning);
+  printf("Calibrating jumping parameters - %d draws\n",tuning);
   begin_time=(int)time((time_t*)NULL);
   AdaptiveMetropolisScale(model,tuning,1000,1,(FILE*)NULL);      // tuning iterations - 1000 iterations before updating - verbose
   end_time=(int)time((time_t*)NULL);
-  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 
   // Reset parametrers
   if (!ReadTransitionMatrices((FILE*)NULL,cmd->parameters_filename_actual,cmd->parameters_header_actual,model) 
       || !Read_VAR_Parameters((FILE*)NULL,cmd->parameters_filename_actual,cmd->parameters_header_actual,model))
-    fprintf(stdout,"Unable to reset parameters after tuning\n");
+    printf("Unable to reset parameters after tuning\n");
 
   // Burn-in period
-  fprintf(stdout,"Burn-in period - %d draws\n",burn_in);
+  printf("Burn-in period - %d draws\n",burn_in);
   for (check=period, count=1; count <= burn_in; count++)
     {
       DrawAll(model);
@@ -109,15 +109,15 @@ int main(int nargs, char **args)
       if (count == check)
 	{
 	  check+=period;
-	  fprintf(stdout,"%d iterations completed out of %d\n",count,burn_in);
+	  printf("%d iterations completed out of %d\n",count,burn_in);
 	}
     }
   end_time=(int)time((time_t*)NULL);
-  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
   ResetMetropolisInformation(p);
  
   // Simulation
-  fprintf(stdout,"Simulating - %d draws\n",iterations);
+  printf("Simulating - %d draws\n",iterations);
   for (check=period, output=thinning, count=1; count <= iterations; count++)
     {
       DrawAll(model);
@@ -136,11 +136,11 @@ int main(int nargs, char **args)
       if (count == check)
 	{
 	  check+=period;
-	  fprintf(stdout,"%d(%d) iterations completed out of %d(%d)\n",count,thinning,iterations,thinning);
+	  printf("%d(%d) iterations completed out of %d(%d)\n",count,thinning,iterations,thinning);
 	}
     }
   end_time=(int)time((time_t*)NULL);
-  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 
   // clean up
   fclose(f_out);
diff --git a/matlab/swz/c-code/sbvar/var/VARbase.c b/matlab/swz/c-code/sbvar/var/VARbase.c
index 11c004d9b8abc68c0a0dd973ca3ebd7b8860c052..f5bf6a4ffbcf7b34f275e869a451eee546200b90 100644
--- a/matlab/swz/c-code/sbvar/var/VARbase.c
+++ b/matlab/swz/c-code/sbvar/var/VARbase.c
@@ -636,7 +636,7 @@ int **ExpandTranslationTable(int **table, TMarkovStateVariable *sv, TMarkovState
   free(idx);
 
   // verbose
-  //dw_PrintArray(stdout,rtable,(char*)NULL); fprintf(stdout,"\n"); dw_PrintArray(stdout,table,(char*)NULL); getchar();
+  //dw_PrintArray(stdout,rtable,(char*)NULL); printf("\n"); dw_PrintArray(stdout,table,(char*)NULL); getchar();
   //dw_PrintArray(stdout,sv->Index,(char*)NULL); getchar();
 
   return rtable;
@@ -1479,7 +1479,7 @@ static void GetProposedJump_A0(TVector b, int j, int k, T_VAR_Parameters *p)
       FreeMatrix(M0);
     }
   S=MatrixInnerProductSymmetric((TMatrix)NULL,p->U[j],YY);
-  //dw_PrintMatrix(stdout,S,"%.17le "); fprintf(stdout,"\n");
+  //dw_PrintMatrix(stdout,S,"%.17le "); printf("\n");
   AddMM(S,S,p->inverse_b0_prior[j]);
 
   //dw_PrintMatrix(stdout,S,"%.17le "); fgetc(stdin);
@@ -1489,20 +1489,20 @@ static void GetProposedJump_A0(TVector b, int j, int k, T_VAR_Parameters *p)
   dw_NormalVector(b);
   if (!InverseProductUV(b,CholeskyUT(S,S),b))
     {
-      fprintf(stdout,"Error in GetProposedJump_A0()\n");
-      fprintf(stdout,"j = %d, k = %d\n,Prior =\n",j,k);
+      printf("Error in GetProposedJump_A0()\n");
+      printf("j = %d, k = %d\n,Prior =\n",j,k);
       dw_PrintMatrix(stdout,p->inverse_b0_prior[j],"%lg ");
-      fprintf(stdout,"S =\n");
+      printf("S =\n");
       dw_PrintMatrix(stdout,S,"%lg ");
       exit(1);
     }
   dw_SetTerminalErrors(terminal_errors);
 /*   else */
 /*     { */
-/*       fprintf(stdout,"GetProposedJump_A0()\n"); */
-/*       fprintf(stdout,"j = %d, k = %d\n,Prior =\n",j,k); */
+/*       printf("GetProposedJump_A0()\n"); */
+/*       printf("j = %d, k = %d\n,Prior =\n",j,k); */
 /*       dw_PrintMatrix(stdout,p->inverse_b0_prior[j],"%lg "); */
-/*       fprintf(stdout,"S =\n"); */
+/*       printf("S =\n"); */
 /*       dw_PrintMatrix(stdout,S,"%lg "); */
 /*       getchar(); */
 /*     } */
@@ -2730,7 +2730,7 @@ void SetLogPriorConstant_VAR(T_VAR_Parameters *p)
       p->log_prior_constant+=p->nvars*log(2);
       break;
     default:
-      fprintf(stdout,"Unknown normalization type\n");
+      printf("Unknown normalization type\n");
       exit(1);
     }
 }
@@ -3037,7 +3037,7 @@ int Normalize_VAR(T_VAR_Parameters *p)
     case VAR_NORMALIZATION_WZ: return WZ_Normalize(p);
     case VAR_NORMALIZATION_NONE: return 0;
     default:
-      fprintf(stdout,"Unknown normalization type\n");
+      printf("Unknown normalization type\n");
       exit(1);
     }
 }
diff --git a/matlab/swz/c-code/sbvar/var/VARio_matlab.c b/matlab/swz/c-code/sbvar/var/VARio_matlab.c
index f4f2a427812b2b7401142daf3384cee259fd0b32..9f8051bb6cfd9bf8999bf3859555506b4742b667 100644
--- a/matlab/swz/c-code/sbvar/var/VARio_matlab.c
+++ b/matlab/swz/c-code/sbvar/var/VARio_matlab.c
@@ -430,13 +430,13 @@ TStateModel* CreateStateModel_VAR_matlab(char *filename)
 /*       YY=TransposeProductMM((TMatrix)NULL,Ui[j],Sigma); */
 /*       ZZ=ProductMM((TMatrix)NULL,YY,Ui[j]); */
 
-/*       fprintf(stdout,"Computed[%d]\n",j); dw_PrintMatrix(stdout,ZZ,"%le "); fprintf(stdout,"\n"); */
-/*       fprintf(stdout,"File[%d]\n",j); dw_PrintMatrix(stdout,H0[j],"%le "); fprintf(stdout,"\n"); */
+/*       printf("Computed[%d]\n",j); dw_PrintMatrix(stdout,ZZ,"%le "); printf("\n"); */
+/*       printf("File[%d]\n",j); dw_PrintMatrix(stdout,H0[j],"%le "); printf("\n"); */
 /*       max=0.0; */
 /*       for (ii=0; ii < RowM(ZZ); ii++) */
 /* 	  for (jj=0; jj < ColM(ZZ); jj++) */
 /* 	    if (max < fabs(ElementM(H0[j],ii,jj) - ElementM(ZZ,ii,jj))) max=fabs(ElementM(H0[j],ii,jj) - ElementM(ZZ,ii,jj)); */
-/*       fprintf(stdout,"H0: max[%d] = %le\n",j,max); */
+/*       printf("H0: max[%d] = %le\n",j,max); */
 
 /*       FreeMatrix(ZZ); */
 /*       FreeMatrix(YY); */
@@ -476,7 +476,7 @@ TStateModel* CreateStateModel_VAR_matlab(char *filename)
 /*       for (ii=0; ii < RowM(ZZ); ii++) */
 /* 	  for (jj=0; jj < ColM(ZZ); jj++) */
 /* 	    if (max < fabs(ElementM(H0[j],ii,jj) - ElementM(ZZ,ii,jj))) max=fabs(ElementM(H0[j],ii,jj) - ElementM(ZZ,ii,jj)); */
-/*       fprintf(stdout,"max[%d] = %le\n",j,max); */
+/*       printf("max[%d] = %le\n",j,max); */
 
 /*       FreeMatrix(ZZ); */
 /*       FreeMatrix(YY); */
diff --git a/matlab/swz/c-code/sbvar/var/estimate.c b/matlab/swz/c-code/sbvar/var/estimate.c
index ab493662741b205ab736a65c9978eec2926ed919..47da3b2b3076b0efa9bf1efea3ca8eea7f1fa844 100644
--- a/matlab/swz/c-code/sbvar/var/estimate.c
+++ b/matlab/swz/c-code/sbvar/var/estimate.c
@@ -289,7 +289,7 @@ int GetLastIteration(FILE *f_in, TStateModel *model, TEstimateInfo *estimate)
 	      sprintf(header=(char*)malloc(strlen(fmt) + i - 1),fmt,k);
 	      if (ReadTransitionMatrices(f_in,(char*)NULL,header,model) && Read_VAR_Parameters(f_in,(char*)NULL,header,model))
 		{
-		  fprintf(stdout,"Using intermediate output - %s\n",header);
+		  printf("Using intermediate output - %s\n",header);
 		  estimate->initialization_header=header;
 		  dw_SetTerminalErrors(terminal_errors);
 		  return 1;
@@ -375,7 +375,7 @@ TStateModel* GetModelFromCommandLine(int nargs, char **args, TEstimateInfo *esti
 	  ReadTransitionMatrices(f_in,(char*)NULL,header,model);
 	  Read_VAR_Parameters(f_in,(char*)NULL,header,model);
 	  fclose(f_in);
-	  fprintf(stdout,"Using final output\n");
+	  printf("Using final output\n");
 	  estimate->specification_filename=filename;
 	  estimate->initialization_filename=filename;
 	  estimate->initialization_header=header;	  
@@ -415,7 +415,7 @@ TStateModel* GetModelFromCommandLine(int nargs, char **args, TEstimateInfo *esti
 	      ReadTransitionMatrices(f_in,(char*)NULL,header,model);
 	      Read_VAR_Parameters(f_in,(char*)NULL,header,model);
 	      fclose(f_in);
-	      fprintf(stdout,"Using initial data\n");
+	      printf("Using initial data\n");
 	      estimate->initialization_filename=filename;
 	      estimate->initialization_header=header;
 	      if (d2) free(d2);	    
@@ -614,7 +614,7 @@ int main(int nargs, char **args)
   //=== Help Screen ===
   if (dw_FindArgument_String(nargs,args,"h") != -1)
     {
-      fprintf(stdout,"print_draws <options>\n");
+      printf("print_draws <options>\n");
       PrintHelpMessages(stdout,include_help,additional_help);
       return 0;
     }
@@ -623,11 +623,11 @@ int main(int nargs, char **args)
   seed=dw_ParseInteger_String(nargs,args,"gs",0);
   dw_initialize_generator(seed);
 
-  fprintf(stdout,"Reading initial data...\n");
+  printf("Reading initial data...\n");
   if (model=SetupFromCommandLine(nargs,args,&estimate))
     {
       // Estimation
-      fprintf(stdout,"Beginning estimation...\n");
+      printf("Beginning estimation...\n");
       begin_time=time((time_t*)NULL);
       FindMode_VAR_csminwel(model,estimate);
       end_time=time((time_t*)NULL);
diff --git a/matlab/swz/c-code/sbvar/var/mhm_VAR.c b/matlab/swz/c-code/sbvar/var/mhm_VAR.c
index 90cf6421a474112a16ffe450a7ec4f0a6ad0cd7e..db2133babfa096e2d95d28f2f7a435d30b62529e 100644
--- a/matlab/swz/c-code/sbvar/var/mhm_VAR.c
+++ b/matlab/swz/c-code/sbvar/var/mhm_VAR.c
@@ -173,7 +173,7 @@ void BurnIn(T_MHM *mhm, int iterations, int period)
 	    }
 
 	  printf("Total Elapsed Time: %d seconds\n",(int)time((time_t*)NULL) - begin_time);
-	  fprintf(stdout,"%d iterations completed out of %d\n",count,iterations);
+	  printf("%d iterations completed out of %d\n",count,iterations);
 	  PrintJumps(stdout,(T_VAR_Parameters*)(mhm->model->theta));
 	}
     }
@@ -204,7 +204,7 @@ void ComputeMeanVariance_MHM(T_MHM *mhm, int iterations, int period)
   InitializeMatrix(mhm->variance,0.0);
 
   // loop and accumulate 1st and 2nd non-central moments
-  fprintf(stdout,"Beginning mean and variance estimation -- %d iterations.\n",iterations);
+  printf("Beginning mean and variance estimation -- %d iterations.\n",iterations);
   begin_time=time((time_t*)NULL);
   for (count=1; count <= iterations; count++)
     {
@@ -227,8 +227,8 @@ void ComputeMeanVariance_MHM(T_MHM *mhm, int iterations, int period)
       if (count == check)
 	{
 	  check+=period;
-	  fprintf(stdout,"Total Elapsed Time: %d seconds\n",(int)time((time_t*)NULL) - begin_time);
-	  fprintf(stdout,"%d iterations completed out of %d\n",count,iterations);
+	  printf("Total Elapsed Time: %d seconds\n",(int)time((time_t*)NULL) - begin_time);
+	  printf("%d iterations completed out of %d\n",count,iterations);
 	  PrintJumps(stdout,(T_VAR_Parameters*)(mhm->model->theta));
 	}
     }
diff --git a/matlab/swz/c-code/sbvar/var/mhm_VAR_main_1.c b/matlab/swz/c-code/sbvar/var/mhm_VAR_main_1.c
index ba0d3eb0d273fea67d079b6ca42c3679d0e6b2ac..2ebc7fce4eba0bddf8c743c7d005b3f506f549f4 100644
--- a/matlab/swz/c-code/sbvar/var/mhm_VAR_main_1.c
+++ b/matlab/swz/c-code/sbvar/var/mhm_VAR_main_1.c
@@ -417,7 +417,7 @@ int main(int nargs, char **args)
 	  initial_time=begin_time=time((time_t*)NULL);
 	  BurnIn_AdaptiveMetropolisScale(mhm,mhm->n_burn1,1000);
 	  end_time=time((time_t*)NULL);
-	  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+	  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 	}
 
       // After first burn-in
@@ -438,10 +438,10 @@ int main(int nargs, char **args)
 	  initial_time=begin_time=time((time_t*)NULL);
 	  BurnIn(mhm,mhm->n_burn2,1000);
 	  end_time=time((time_t*)NULL);
-	  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+	  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 
-	  fprintf(stdout,"Number inconsistent normalizations: %d\n",((T_VAR_Parameters*)mhm->model->theta)->WZ_inconsistancies);
-	  fprintf(stdout,"Number singular inverse variances: %d\n\n",Get_VAR_Improper_Distribution_Counter());
+	  printf("Number inconsistent normalizations: %d\n",((T_VAR_Parameters*)mhm->model->theta)->WZ_inconsistancies);
+	  printf("Number singular inverse variances: %d\n\n",Get_VAR_Improper_Distribution_Counter());
 	}
 
       fclose(f_out_intermediate_draws);
@@ -463,20 +463,20 @@ int main(int nargs, char **args)
 	  begin_time=time((time_t*)NULL);
 	  ComputeMeanVariance_MHM(mhm,mhm->n_mean_variance,10000);
 	  end_time=time((time_t*)NULL);
-	  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+	  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 
-	  fprintf(stdout,"Number inconsistent normalizations: %d\n",((T_VAR_Parameters*)mhm->model->theta)->WZ_inconsistancies);
-	  fprintf(stdout,"Number singular inverse variances: %d\n\n",Get_VAR_Improper_Distribution_Counter());
+	  printf("Number inconsistent normalizations: %d\n",((T_VAR_Parameters*)mhm->model->theta)->WZ_inconsistancies);
+	  printf("Number singular inverse variances: %d\n\n",Get_VAR_Improper_Distribution_Counter());
 	}
 
       // Set center to mean if necessary
       if (dw_FindArgument_String(nargs,args,"cm") >= 0)
 	{
-	  fprintf(stdout,"Using mean for center\n");
+	  printf("Using mean for center\n");
 	  mhm->center=mhm->mean;
 	}
       else
-	fprintf(stdout,"Using posterior mode for center\n");
+	printf("Using posterior mode for center\n");
 
 
       // After mean-variance estimation
@@ -507,9 +507,9 @@ int main(int nargs, char **args)
       begin_time=time((time_t*)NULL);
       ComputeModifiedHarmonicMean(mhm,10000);
       end_time=time((time_t*)NULL);
-      fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
-      fprintf(stdout,"Number inconsistent normalizations: %d\n",((T_VAR_Parameters*)mhm->model->theta)->WZ_inconsistancies);
-      fprintf(stdout,"Number singular inverse variances: %d\n\n",Get_VAR_Improper_Distribution_Counter());
+      printf("Elapsed Time: %d seconds\n",end_time - begin_time);
+      printf("Number inconsistent normalizations: %d\n",((T_VAR_Parameters*)mhm->model->theta)->WZ_inconsistancies);
+      printf("Number singular inverse variances: %d\n\n",Get_VAR_Improper_Distribution_Counter());
 
       fclose(mhm->f_out);
 
diff --git a/matlab/swz/c-code/sbvar/var/mhm_VAR_main_2.c b/matlab/swz/c-code/sbvar/var/mhm_VAR_main_2.c
index eeb423b2aa3a0eb10fd26ca7e89c2b16f443c6a2..4169ff3a1e4b817a06816a439850be2ff25b7e4e 100644
--- a/matlab/swz/c-code/sbvar/var/mhm_VAR_main_2.c
+++ b/matlab/swz/c-code/sbvar/var/mhm_VAR_main_2.c
@@ -696,7 +696,7 @@ TVector Create_q(int ndraws, TVector center, TMatrix scale, TStateModel* model,
 
   // timings
   end_time=time((time_t*)NULL);
-  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 
   FreeVector(free_parameters);
 
@@ -779,7 +779,7 @@ TMatrix CreateProposal(int ndraws, TVector center, TMatrix scale, TStateModel* m
 
   // timings
   end_time=time((time_t*)NULL);
-  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 
   FreeVector(free_parameters);
 
@@ -859,7 +859,7 @@ TMatrix CreateProposal_Radius(int ndraws, TVector center, TMatrix scale, TStateM
 
   // timings
   end_time=time((time_t*)NULL);
-  fprintf(stdout,"Elapsed Time: %d seconds\n",end_time - begin_time);
+  printf("Elapsed Time: %d seconds\n",end_time - begin_time);
 
   FreeVector(free_parameters);
 
diff --git a/matlab/swz/c-code/utilities/TZCcode/kalman.c b/matlab/swz/c-code/utilities/TZCcode/kalman.c
index c9804690f6d3152e976696820829a49f3f09ca58..e20a2c192dbf2d435057419ca68848f732756b41 100644
--- a/matlab/swz/c-code/utilities/TZCcode/kalman.c
+++ b/matlab/swz/c-code/utilities/TZCcode/kalman.c
@@ -777,7 +777,7 @@ double tz_kalfiltv(struct TSkalfiltv_tag *kalfiltv_ps)
       }
       else
       {
-         fprintf(stdout, "Fatal error: tz_kalfiltv() in kalman.c: the system is non-stationary solutions\n"
+         printf("Fatal error: tz_kalfiltv() in kalman.c: the system is non-stationary solutions\n"
                          "    and the initial conditions must be supplied by, say, input arguments");
          fflush(stdout);
          exit( EXIT_FAILURE );
@@ -1311,9 +1311,9 @@ int InitializeKalman_z10_P10(struct TSkalfilmsinputs_1stapp_tag *kalfilmsinputs_
                      fprintf(FPTR_DEBUG, ".../kalman.c/InitializeKalman_z10_P10(): the system is non-stationary solutions\n"
                                          "    and see p.003 in LWZ Model II");
                      #else
-                     fprintf(stdout, "\n-----------------\n");
-                     fprintf(stdout, "\nIn grand regime sti=%d\n", sti);
-                     fprintf(stdout, ".../kalman.c/InitializeKalman_z10_P10(): the system is non-stationary solutions\n"
+                     printf("\n-----------------\n");
+                     printf("\nIn grand regime sti=%d\n", sti);
+                     printf(".../kalman.c/InitializeKalman_z10_P10(): the system is non-stationary solutions\n"
                                      "    and see p.003 in LWZ Model II");
                      #endif
                   }
@@ -1836,8 +1836,8 @@ double tz_logTimetCondLH_kalfilms_1st_approx(int st, int inpt, struct TSkalfilms
             }
             else
             {
-               fprintf(stdout, "\n-----------------\n");
-               fprintf(stdout, "\nIn regime comsti_c=%d and sti_v=%d and at time=%d\n", comsti_c, sti_v, 0);
+               printf("\n-----------------\n");
+               printf("\nIn regime comsti_c=%d and sti_v=%d and at time=%d\n", comsti_c, sti_v, 0);
                fn_DisplayError("kalman.c/tz_Refresh_z_T7P_T_in_kalfilms_1st_approx(): the system is non-stationary solutions\n"
                              "    and the initial conditions must be supplied by, say, input arguments");
                fflush(stdout);
@@ -2464,8 +2464,8 @@ void tz_Refresh_z_T7P_T_in_kalfilms_1st_approx(struct TSkalfilmsinputs_tag *kalf
          }
          else
          {
-            fprintf(stdout, "\n-----------------\n");
-            fprintf(stdout, "\nIn regime sti_c=%d and sti_v=%d and at time=%d\n", sti_c, sti_v, 0);
+            printf("\n-----------------\n");
+            printf("\nIn regime sti_c=%d and sti_v=%d and at time=%d\n", sti_c, sti_v, 0);
             fn_DisplayError("kalman.c/tz_Refresh_z_T7P_T_in_kalfilms_1st_approx(): the system is non-stationary solutions\n"
                           "    and the initial conditions must be supplied by, say, input arguments");
             fflush(stdout);
diff --git a/matlab/swz/c-code/utilities/TZCcode/mathlib.c b/matlab/swz/c-code/utilities/TZCcode/mathlib.c
index 9847519bf3ce66311a1fb76848c051726ee28722..b2728b82a3907ac5d99a0437b13335758ff0d128 100644
--- a/matlab/swz/c-code/utilities/TZCcode/mathlib.c
+++ b/matlab/swz/c-code/utilities/TZCcode/mathlib.c
@@ -4304,7 +4304,7 @@ void TransposeSquare(TSdmatrix *B_dm, TSdmatrix *A_dm) {
       fprintf(FPTR_DEBUG, "\nWARNING: .../mathlib.c/TransposeSquare():  the matrix is already both SU and SL, so there is no need to transpose.\n");
       fflush(FPTR_DEBUG);
       #else
-      fprintf(stdout, "\nWARNING: .../mathlib.c/TransposeSquare():  the matrix is already both SU and SL, so there is no need to transpose.\n");
+      printf("\nWARNING: .../mathlib.c/TransposeSquare():  the matrix is already both SU and SL, so there is no need to transpose.\n");
       fflush(stdout);
       #endif
 
diff --git a/matlab/swz/c-code/utilities/TZCcode/tzmatlab.c b/matlab/swz/c-code/utilities/TZCcode/tzmatlab.c
index c27b731f8b2076a51bd4ce20a8e53955f585c1ba..3876c3fa94e57b9dadab88a9386c0bf2d20935c4 100644
--- a/matlab/swz/c-code/utilities/TZCcode/tzmatlab.c
+++ b/matlab/swz/c-code/utilities/TZCcode/tzmatlab.c
@@ -3,7 +3,7 @@
    fprintf(FPTR_DEBUG, "\nWARNING: .../mathlib.c/TransposeSquare():  the matrix is already both SU and SL, so there is no need to transpose.\n");
    fflush(FPTR_DEBUG);
    #else
-   fprintf(stdout, "\nWARNING: .../mathlib.c/TransposeSquare():  the matrix is already both SU and SL, so there is no need to transpose.\n");
+   printf("\nWARNING: .../mathlib.c/TransposeSquare():  the matrix is already both SU and SL, so there is no need to transpose.\n");
    fflush(stdout);
    #endif
 /**/
@@ -48,7 +48,7 @@ void fn_DisplayError(char *msg_s)
              "  %s!\n", msg_s);
       fflush(FPTR_DEBUG);
    #else
-      fprintf(stdout, "\nFatal Error:\n"
+      printf("\nFatal Error:\n"
              "\t %s!\n", msg_s);
       fflush(stdout);
    #endif