diff --git a/prism/src/prism/ProbModelChecker.java b/prism/src/prism/ProbModelChecker.java
index 2cc5221a..e820e94d 100644
--- a/prism/src/prism/ProbModelChecker.java
+++ b/prism/src/prism/ProbModelChecker.java
@@ -45,6 +45,7 @@ import jdd.JDDNode;
 import jdd.JDDVars;
 import mtbdd.PrismMTBDD;
 import parser.ast.Expression;
+import parser.ast.ExpressionFunc;
 import parser.ast.ExpressionProb;
 import parser.ast.ExpressionReward;
 import parser.ast.ExpressionSS;
@@ -861,7 +862,7 @@ public class ProbModelChecker extends NonProbModelChecker
 			return checkRewardReach((ExpressionTemporal) expr, stateRewards, transRewards);
 		}
 		else if (Expression.isCoSafeLTLSyntactic(expr)) {
-			throw new PrismNotSupportedException("Co-safe reward properties not yet supported for DTMCs in this engine");
+			return checkRewardCoSafeLTL(expr, stateRewards, transRewards);
 		}
 		throw new PrismException("R operator contains a path formula that is not syntactically co-safe: " + expr);
 	}
@@ -899,6 +900,151 @@ public class ProbModelChecker extends NonProbModelChecker
 		return rewards;
 	}
 
+	// co-safe LTL reward
+
+	protected StateValues checkRewardCoSafeLTL(Expression expr, JDDNode stateRewards, JDDNode transRewards) throws PrismException
+	{
+		LTLModelChecker mcLtl;
+		StateValues rewardsProduct = null, rewards = null;
+		Expression ltl;
+		Vector<JDDNode> labelDDs;
+		DA<BitSet, ? extends AcceptanceOmega> da;
+		ProbModel modelProduct;
+		ProbModelChecker mcProduct;
+		JDDNode startMask;
+		JDDVars daDDRowVars, daDDColVars;
+		int i;
+		long l;
+
+		if (Expression.containsTemporalTimeBounds(expr)) {
+			if (model.getModelType().continuousTime()) {
+				throw new PrismException("DA construction for time-bounded operators not supported for " + model.getModelType()+".");
+			}
+
+			if (expr.isSimplePathFormula()) {
+				// Convert simple path formula to canonical form,
+				// DA is then generated by LTL2RabinLibrary.
+				//
+				// The conversion to canonical form has to happen here, because once
+				// checkMaximalStateFormulas has been called, the formula should not be modified
+				// anymore, as converters may expect that the generated labels for maximal state
+				// formulas only appear positively
+				expr = Expression.convertSimplePathFormulaToCanonicalForm(expr);
+			} else {
+				throw new PrismException("Time-bounded operators not supported in LTL: " + expr);
+			}
+		}
+
+		// Can't do "dfa" properties yet
+		if (expr instanceof ExpressionFunc && ((ExpressionFunc) expr).getName().equals("dfa")) {
+			throw new PrismException("Model checking for \"dfa\" specifications not supported yet");
+		}
+
+		// For LTL model checking routines
+		mcLtl = new LTLModelChecker(prism);
+
+		// Model check maximal state formulas
+		labelDDs = new Vector<JDDNode>();
+		ltl = mcLtl.checkMaximalStateFormulas(this, model, expr.deepCopy(), labelDDs);
+
+		// Convert LTL formula to deterministic automaton (DA)
+		mainLog.println("\nBuilding deterministic automaton (for " + ltl + ")...");
+		l = System.currentTimeMillis();
+		LTL2DA ltl2da = new LTL2DA(prism);
+		AcceptanceType[] allowedAcceptance = {
+				AcceptanceType.RABIN,
+				AcceptanceType.REACH
+		};
+		da = ltl2da.convertLTLFormulaToDA(ltl, constantValues, allowedAcceptance);
+		mainLog.println(da.getAutomataType()+" has " + da.size() + " states, " + da.getAcceptance().getSizeStatistics()+".");
+		l = System.currentTimeMillis() - l;
+		mainLog.println("Time for deterministic automaton translation: " + l / 1000.0 + " seconds.");
+		// If required, export DA 
+		if (prism.getSettings().getExportPropAut()) {
+			mainLog.println("Exporting DA to file \"" + prism.getSettings().getExportPropAutFilename() + "\"...");
+			PrismLog out = new PrismFileLog(prism.getSettings().getExportPropAutFilename());
+			da.print(out, prism.getSettings().getExportPropAutType());
+			out.close();
+			//da.printDot(new java.io.PrintStream("da.dot"));
+		}
+
+		// Build product of Markov chain and automaton
+		// (note: might be a CTMC - StochModelChecker extends this class)
+		mainLog.println("\nConstructing MC-"+da.getAutomataType()+" product...");
+		daDDRowVars = new JDDVars();
+		daDDColVars = new JDDVars();
+		l = System.currentTimeMillis();
+		modelProduct = mcLtl.constructProductMC(da, model, labelDDs, daDDRowVars, daDDColVars);
+		l = System.currentTimeMillis() - l;
+		mainLog.println("Time for product construction: " + l / 1000.0 + " seconds.");
+		mainLog.println();
+		modelProduct.printTransInfo(mainLog, prism.getExtraDDInfo());
+		// Output product, if required
+		if (prism.getExportProductTrans()) {
+			try {
+				mainLog.println("\nExporting product transition matrix to file \"" + prism.getExportProductTransFilename() + "\"...");
+				modelProduct.exportToFile(Prism.EXPORT_PLAIN, true, new File(prism.getExportProductTransFilename()));
+			} catch (FileNotFoundException e) {
+				mainLog.printWarning("Could not export product transition matrix to file \"" + prism.getExportProductTransFilename() + "\"");
+			}
+		}
+		if (prism.getExportProductStates()) {
+			mainLog.println("\nExporting product state space to file \"" + prism.getExportProductStatesFilename() + "\"...");
+			PrismFileLog out = new PrismFileLog(prism.getExportProductStatesFilename());
+			modelProduct.exportStates(Prism.EXPORT_PLAIN, out);
+			out.close();
+		}
+
+		// Adapt reward info to product model
+		JDD.Ref(stateRewards);
+		JDD.Ref(modelProduct.getReach());
+		JDDNode stateRewardsProduct = JDD.Apply(JDD.TIMES, stateRewards, modelProduct.getReach());
+		JDD.Ref(transRewards);
+		JDD.Ref(modelProduct.getTrans01());
+		JDDNode transRewardsProduct = JDD.Apply(JDD.TIMES, transRewards, modelProduct.getTrans01());
+		
+		// Find accepting states + compute reachability rewards
+		AcceptanceOmegaDD acceptance = da.getAcceptance().toAcceptanceDD(daDDRowVars);
+		JDDNode acc = null;
+		if (acceptance instanceof AcceptanceReachDD) {
+			mainLog.println("\nSkipping BSCC computation since acceptance is defined via goal states...");
+			acc = ((AcceptanceReachDD) acceptance).getGoalStates();
+			JDD.Ref(modelProduct.getReach());
+			acc = JDD.And(acc, modelProduct.getReach());
+		} else {
+			mainLog.println("\nFinding accepting BSCCs...");
+			acc = mcLtl.findAcceptingBSCCs(acceptance, modelProduct);
+		}
+		acceptance.clear();
+		mainLog.println("\nComputing reachability rewards...");
+		mcProduct = createNewModelChecker(prism, modelProduct, null);
+		rewardsProduct = mcProduct.computeReachRewards(modelProduct.getTrans(), modelProduct.getTrans01(), stateRewardsProduct, transRewardsProduct, acc);
+
+		// Convert reward vector to original model
+		// First, filter over DRA start states
+		startMask = mcLtl.buildStartMask(da, labelDDs, daDDRowVars);
+		JDD.Ref(model.getReach());
+		startMask = JDD.And(model.getReach(), startMask);
+		rewardsProduct.filter(startMask);
+		// Then sum over DD vars for the DRA state
+		rewards = rewardsProduct.sumOverDDVars(daDDRowVars, model);
+
+		// Deref, clean up
+		JDD.Deref(stateRewardsProduct);
+		JDD.Deref(transRewardsProduct);
+		rewardsProduct.clear();
+		modelProduct.clear();
+		for (i = 0; i < labelDDs.size(); i++) {
+			JDD.Deref(labelDDs.get(i));
+		}
+		JDD.Deref(acc);
+		JDD.Deref(startMask);
+		daDDRowVars.derefAll();
+		daDDColVars.derefAll();
+
+		return rewards;
+	}
+
 	// steady state reward
 
 	protected StateValues checkRewardSS(ExpressionTemporal expr, JDDNode stateRewards, JDDNode transRewards) throws PrismException